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Abstract 

We present a quantitative continuum theory of "flocking": the collective 
coherent motion of large numbers of self-propelled organisms. In agreement 
with everyday experience, our model predicts the existence of an "ordered 
phase" of flocks, in which all members of even an arbitrarily large flock move 
together with the same mean velocity (w) 7^ 0. This coherent motion of the 
flock is an example of spontaneously broken symmetry: no preferred direction 
for the motion is picked out a priori in the model; rather, each flock is allowed 
to, and does, spontaneously pick out some completely arbitrary direction to 
move in. By analyzing our model we can make detailed, quantitative pre- 
dictions for the long-distance, long-time behavior of this "broken symmetry 
state". The "Goldstone modes" associated with this "spontaneously broken 
rotational symmetry" are fluctuations in the direction of motion of a large 
part of the flock away from the mean direction of motion of the flock as a 
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whole. These "Goldstone modes" mix with modes associated with conser- 
vation of bird number to produce propagating sound modes. These sound 
modes lead to enormous fluctuations of the density of the flock, far larger, at 
long wavelengths, than those in, e. g., an equilibrium gas. Our model is sim- 
ilar in many ways to the Navier-Stokes equations for a simple compressible 
fluid; in other ways, it resembles a relaxational time dependent Ginsburg- 
Landau theory for an n = d component isotropic ferromagnet. In spatial 
dimensions d > 4, the long distance behavior is correctly described by a 
linearized theory, and is equivalent to that of an unusual but nonetheless 
equilibrium model for spin systems. For d < A, non-linear fluctuation effects 
radically alter the long distance behavior, making it different from that of 
any known equilibrium model. In particular, we find that in d = 2, where 
we can calculate the scaling exponents exactly, flocks exhibit a true, long- 
range ordered, spontaneously broken symmetry state, in contrast to equilib- 
rium systems, which cannot spontaneously break a continuous symmetry in 
d = 2 (the "Mermin- Wagner" theorem). We make detailed predictions for 
various correlation functions that could be measured either in simulations, or 
by quantitative imaging of real flocks. We also consider an anisotropic model, 
in which the birds move preferentially in an "easy" (e. g., horizontal) plane, 
and make analogous, but quantitatively different, predictions for that model 
as well. For this anisotropic model, we obtain exact scaling exponents for all 
spatial dimensions, including the physically relevant case d = 3. 
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I. INTRODUCTION 



A wide variety of non-equilibrium dynamical systems with many degrees of freedom 
have recently been studied using powerful techniques developed for equilibrium condensed 
matter physics (e.g., scaling, the renormalization group, etc). One of the most familiar 
examples of a many degree of freedom, non-equilibrium dynamical system is a large flock 
of birds. Myriad other examples of the collective, coherent motion of large numbers of self- 
propelled organisms occur in biology: schools of fish, swarms of insects, slime molds, herds 
of wildebeest, etc. 

Recently, a number of simulations of this phenomenon have been performed. Fol- 
lowing Reynolds [§, we will use the term "bold" and bird interchangeably for the particles 
in these simulations. All of these simulations have several essential features in common: 

1. A large number (a "flock") of point particles ("boids") each move over time through a 
space of dimension d (= 2,3,...), attempting at all times to "follow" (i.e., move in the 
same direction as) its neighbors. 

2. The interactions are purely short ranged: each "bold" only responds to its neighbors, 
defined as those "boids" within some fixed, finite distance Ro, which is assumed to be 
much less than L, the size of the "flock." 

3. The "following" is not perfect: the "boids" make errors at all times, which are modeled 
as a stochastic noise. This noise is assumed to have only short ranged spatio-temporal 
correlations. 

4. The underlying model has complete rotational symmetry: the flock is equally likely, a 
priori, to move in any direction. 

The development of a non-zero mean center of mass velocity (v) for the flock as a whole 
therefore requires spontaneous breaking of a continuous symmetry (namely, rotational). 
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In an earlier paper [Q, we formulated a continuum model for such dynamics of flocking, 
and obtained some exact results for that model in spatial dimensions d = 2 (appropriate for 
the description of the motion of land animals on the earth's surface). Our most surprising 
result 0] was that two-dimensional moving herds with strictly short-ranged interactions 
appear to violate the Mermin- Wagner theorem in that they can acquire long-ranged 
order, by picking out a consistent direction of motion across an arbitrarily large herd, despite 
the fact that this involves spontaneously breaking a continuous (rotational) symmetry. 

Of course, this result does not, in fact, violate the Mermin- Wagner theorem, since flocks 
are a non-equilibrium dynamical system. What is fascinating (at least to us) about our result 
is that the non-equilibrium aspects of the flock dynamics that make the long-distance, long- 
time behavior of the flock different from that of otherwise analogous equilibrium systems 
are fundamentally non-linear, strong-fluctuation effects. Indeed, a "breakdown of linearized 
hydrodynamics," analogous to that long known to occur in equilibrium fluids in spatial 
dimensions d = 2, occurs in flocks for all d < A. This breakdown of linearized hydrodynamics 
is essential to the very existence of the ordered state in d = 2. Furthermore, it has dramatic 
consequences even ioi d > 2. 

The physics of this breakdown is very simple: above d = 4, where the breakdown does 
not occur, information about what is going on in one part of the flock can be transmitted to 
another part of the flock only by being passed sequentially through the intervening neighbors 
via the assumed short ranged interactions. Below d = 4, where the breakdown occurs, 
this slow, diffusive transport of information is replaced by direct, convective transport: 
fluctuations in the local velocity of the flock became so large, in these lower dimensions, 
that the motion of one part of the flock relative to another becomes the principle means of 
information transport, because it becomes faster than diffusion. There is a sort of "negative 
feedback," in that this improved transport actually suppresses the very fluctuations that 
give rise to it, leading to long-ranged order in = 2. The purpose of the present paper is to 
study the properties of the "ordered state" of the flock, i.e., the state in which all members 
of the flock are moving in the same average direction. Specifically, we will: 
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1. give the details of the derivation of the results of reference [Q, and give detailed pre- 
dictions for numerous correlation functions that can be measured in both experiments 
and simulations. In particular, we will show that two propagating sound modes exist 
in flocks, with unusually anisotropic speeds, whose detailed dependence on the direc- 
tion of propagation we predict, making possible extremely stringent quantitative tests 
of our theory. We also calculate their attenuations, which show highly anomalous, and 
strongly anisotropic, scaling, 

2. formulate and study the most complete generalization of the model of reference [|] for 
spatial dimensions d > 2, and 

3. include the effect of spatial anisotropy (e.g., the fact that birds prefer to fly horizontally 
rather than vertically) on flock motion. 

We describe the flock with coarse grained density and velocity fields p{r,t) and v{r,t), 
respectively, giving the average number density and velocity of the birds at time t within 
some coarse graining distance io of a given position r in space. The coarse graining distance 
io is chosen to be as small as possible, consistent with being large enough that the averaging 
can be done sensibly (in particular, must be greater than the mean interbird distance). 
Our description is then valid for distances large compared to io, and for times t much greater 
than some microscopic time to, presumably of order io/vr, where vt is a typical speed of 
a bird. Collective motion of the flock as a whole then requires that {v{f,t)) ^ 0; where 
the averaging can be considered an ensemble average, a time average, or a spatial average. 
Equivalently, long ranged order must develop for the flock as a whole to move; i.e., the 
equal-time velocity auto-correlation function: 





must approach a non-zero constant as the separation \R 



oo; specifically: 
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Thus the average velocity (v) of the flock is precisely analogous to the order parameter 
(s) in a ferromagnetic system, where s is a local spin. 

Our most dramatic result is that an intrinsically non-equilibrium and nonlinear feature of 
our model, namely, convection, suppresses fluctuations of the velocity v at long wavelengths, 
making them much smaller than the analogous s fluctuations found in ferromagnets, for 
all spatial dimensions of the flock d < 4. Specifically, the connected piece Cc{R) of the 
correlation function C{R), defined as: 

Cc{R) ^C(r)- hm C (r') , (1.3) 

— * 

which is a measure of the fluctuations, decays to zero much more rapidly, as |it!| — > oo, than 
the analogous correlation function in magnets. Quantitatively, for points whose separation 
R = R± lies perpendicular to the mean direction of motion of the flock, 

Cc{R) oc (1.4) 

where the universal "roughness exponent" 

X = -l (1-5) 

exactly, in d = 2, and is < 1 — |, its value in magnetic systems, for all d < 4. For d > i, 
X = 1 — I for flocks as well as for magnets. 

The physical mechanism for this suppression of fluctuations is easy to understand: in- 
creased fluctuations in the direction of motion of different parts of the flock actually enhance 
the exchange of information between those different parts. This exchange, in turn, suppresses 
those very fluctuations, since the interactions between birds tend to make them all move in 
the same direction. 

These non-equilibrium effects also lead to a spatial anisotropy of scaling between the 
direction along (||) and those orthogonal to (±) the mean velocity {v). The physical origin 
of the anisotropy is also simple: if birds make small errors 66 in their direction of motion, 
their random motion perpendicular to the mean direction of motion {v) is much larger than 
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that along (v); the former being oc S6, while the later is proportional to 1 — cos 5^ ~ 66"^. 
As a result, any equal-time correlation function in the system of any combination of fields 

crosses over from dependence purely on \R±\ to dependence purely on R\\ when 

C 

(1.6) 



\Ri 



where £o is the bird interaction range. 
The universal anisotropy exponent 
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exactly in d = 2, and is < 1 for all d < 4. 

In particular, the connected, equal-time, velocity autocorrelation function Cc{R) obeys 
the scaling law 

^ (^11 /4) 



Cc{r) 



Ri 



2X 



{\R±\/^o) 



X 



:i.8) 



where fv{x) is a universal scaling function. We have, alas, been unable to calculate this 
scaling function, even in d = 2 where we know the exponents exactly. However, the scaling 
form ( |1.8|) immediately implies that 



Cc i R) OC R, 



when R\\/io > ^Rjio) 



[1.9] 



So far, our discussion has focussed on velocity fluctuations. The density p{R, t) shows 
huge fluctuations as well: indeed, at long wavelengths, the fluctuations of the density of birds 
in a flock become infinitely bigger than those in a fluid or an ideal gas. This fact is obvious 
to the eye in a picture of a flock (see Fig. p. Quantitatively, we predict that the spatially 
Fourier transformed, equal time density-density correlation function Cp{q) = (|p(g, t)p) 



obeys the scaling law: 



3-d-(-2x 



,(g±^o) 



Y{9^) oc 



1- 
(1± 



-C-2X 



-2 J-fl!-C-2x 



(4g±)^»^ogii (1-10) 



^11 
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-3+ 



l-d-2x 



"ql, (g±4)^<gii4 



where Y{9^ is a finite, non- vanishing, function of the angle % between the wavevec- 
tor q and the direction of mean fiock motion, q\\ and are the wavevectors parallel and 
perpendicular to the broken symmetry direction, and q±_ = \q±\. 
In d — 2, ( — ^ and x — ~|) so 



^iT'^I, > ioq\\ > (1-11) 

Qw^qI, {q±hf' < q\\h 

The most important thing to note about Cp{q) is that it diverges as |g| — > 0, unlike Cp{q) 
for, say, a simple fluid or gas, or, indeed, for any equilibrium condensed matter system, which 
goes to a finite constant (the compressibility) as \q\ 0. 

This correlation function should be extremely easy to measure in simulations, and in 
experiments on real herds or fiocks, in which, say, video tape allows one to measure the 
positions fi{t) of all the birds (labeled by i) in the fiock at a variety of times t. The recipe 
is simple: 

1. Calculate the complex numbers 

p(g,t) = ^e^«"'-'W (1.12) 

i 

for a variety of ^s. 

2. Average the squared magnitude of this number over time. The result is Cp{q). 

Time dependent correlation functions of p and v in fiocks also show interesting anoma- 
lous scaling behavior. However, it is not so simple to summarize as the equal time cor- 
relation functions. Indeed, time dependent correlation functions (or, equivalcntly, their 
spatio-temporal Fourier transforms) do not have a simple scaling form. This is because the 
collective normal modes of the fiock consist of propagating, damped longitudinal "sound" 
modes (i.e., density waves), as well as, in > 2, shear modes. The sound modes exhibit two 
different types of scaling: the period T of a wave is proportional to its wavelength A (the 



constant of proportionality being the inverse sound speed); while the lifetime r of the mode 
is proportional to A^, with z being another universal exponent. In most systems (e.g., fluids, 
crystals) exhibiting sound modes, z = 2, corresponding to conventional diffusive or viscous 
damping 0. In flocks, however, we find 
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z = -; d = 2 (1.13) 
5 

and z < 2 for all d < 4, for sound modes propagating orthogonal to the mean direction of 
flock motion. That is, sound modes are much more heavily damped at long wavelengths in 
flocks than in most equilibrium condensed matter systems. 
The full dispersion relation for the sound modes is 



oo± = c^{e^)q - tqlf± j (1.14) 

where % is the angle between q and (v), and the direction-dependent sound speeds c± [9^) 
are given by equation ( [4.11| ) of section (IV), with 7 and ai flock-dependent parameters and 
Po the mean number density of "birds" in the flock. A polar plot of these sound speeds is 
given in Fig. (^. The exponent C, is the universal anisotropy exponent described earlier, 
and f±{x) are universal scaling functions which, alas, we have been unable to calculate. 
However, we do know some of their limits: 

f±{x 0) ^ constant > 0; f±{x 00) oc . (1-15) 

Note that this last result implies that the lifetime r of the wave is oc for q\\(.Q 3> 
(fjx^o)^ and |g| — 0; this only happens for directions of propagation very nearly parallel to 
{v). Note also that in d = 2 where = | and C = |, | = 2 and the damping is conventional 
for sound modes propagating parallel to the mean motion of the flock. For all other directions 
of propagation, however, it is unconventional, and characterized by z = |. 

This behavior of the damping (i.e. /ma;), is summarized in Fig. (^. For d > 2, the 
"hydrodynamic" mode structure also includes d — 2 "hyper diffusive" shear modes, with 
identical dispersion relations 



g||^0 



u;s = iq\\-iqlfs[y^A ■ (1.16) 



The dispersion relations for uj± and Ug can be directly probed by measuring the spatio- 
temporally Fourier transformed density-density and velocity-velocity auto-correlation func- 
tions 



Cp{q,u)^{\p{lu)\^) (1.17) 
Cij{q,uj) = {Vi{q,uj)vj{-q,-uj)) 

respectively. Experimentally, or in simulations, Cp [q, uj) can be calculated by temporally 
Fourier transforming the spatially Fourier transformed density ( 1.12|) : 



(n+l)T 

Pn{q,uj)= E Pn{q,t)e-'^\ n = 0,l,2,... (1.18) 

t=nT 

over a set of long "bins" of time intervals of length r ^ to (the "microscopic" time step), 
and then averaging the squared magnitude \p {q,uj)f over bins: 

(1.19) 

Closed form expressions for these correlation functions in terms of the scaling functions 
/i and fs are given in section V. Although these expressions look quite complicated, the 
behavior they predict is really quite simple, as illustrated in Fig. (|), where Cp is plotted 
as a function of uj for fixed q. As shown there, Cp has two sharp peaks at u; = c± [dg-) q, of 



width oc q^fi {j^eY) height oc qj_ 9 {j^eYj ' '^^ ^^^^ simply 

extracted from the position of the peaks, while the exponents z and ( can be determined 
by comparing their widths and heights for different g 's. 

The scaling properties of the flock are completely summarized by the universal exponents 
z, (, and X- In c? = 2, our predictions for these exponents are: 

z^l C-l X^-i (1.20) 

These results are exact and universal for all flocks with the simple symmetries we discussed 
at the outset. 
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For (i 7^ 2, the situation is less clear. We have performed a one loop, 4 — e expansion to 
attempt to calculate these exponents, and find that, to this order, the model appears to have 
a fixed line with continuously varying exponents and x- Whether this is an artifact 
of our one loop calculation, or actually happens, is unclear. A two loop calculation might 
clarify matters, but would be extremely long and tedious. (One loop was hard enough.) 

The origin of this complication is an additional convective non-linearity P] not discussed 
in Ref. This new term (whose coefficient is a parameter we call A2) is unrenormalized 
at one loop order, leading to the apparent fixed line at that order. In two dimensions, this 
extra term can be written as a total derivative, and can be absorbed into the non-linear 
term considered in Ref. [Q. Hence, in c? = 2, the results of [0] are sound. In d > 2, however, 
this new term has a different structure, and could, if it does not renormalize to zero, change 
the exponents Xi a-^^d Since A2 does not renormalize at one loop order, all we can say 
at this point is that there are three possibilities: 

1. At higher order, A2 renormalizes to zero, //this is the case, we can show that 

C^^. X^^. (1.21) 

5 5 5 

exactly, for all d in the range 2 < d < A. Note that these results linearly interpolate 
between the equilibrium results 2; = 2, ^ = 1, and x=l~fin(i = 4, and our 2d 
results 2; = |, = |, and x = ^| in d = 2. 

2. At higher order, A2 grows upon renormalization and reaches a non-zero fixed point 
value A2 at some new fixed point that differs from the A2 = fixed point we've studied 
previously, at which Eqn. ( 1.21|) holds. The exponents x^ ^1 and C would still be 



universal (i. e., depend only on the dimension of space d) for all fiocks in this case, 
but those universal values would be different from Eqn. ( |1.21 ). 



3. A2 is unrenormalized to all orders. Should this happen, A2 would parameterize a fixed 
/me, with continuously varying values of the exponents x, and C- 
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We reiterate: we do not know which of the above possibihties holds for d > 2. However, 
whichever holds is universal] that is, only one of the three possibilities above applies to all 
flocks. We don't, however, know which one that is. 

We also study an anisotropic model for flocking, which incorporates the possibility that 
birds are averse to flying in certain directions (i.e., straight up or straight down). In partic- 
ular, we consider the case in which, for arbitrary spatial dimensions d > 2, there is an easy 
plane for motion (i.e., a. de = 2 dimensional subspace of the full d-dimensional space). In 
this case, the relevant pieces of the A2 vertex become a total derivative, and can be absorbed 
into the non-linear term considered in reference 0], for all spatial dimensions d, not just 
= 2 as in the isotropic model. Hence, we are able to obtain exact exponents for this 
problem for all spatial dimensions d, not just d = 2. 

We again find anisotropic, anomalous scaling for d < 4. The anisotropy of scaling is 
between the direction in the easy plane (call it x) perpendicular to the mean direction of 
motion (call it y; which, of course, also lies in the easy plane), and all rf — 1 other directions, 
including y (see figure 5). That is, the equal time, velocity- velocity autocorrelation function 
obeys the scaling law: 



[1.22) 



where th denotes the d — 2 components of r in the "hard" directions orthogonal to the easy 
plane, with the scaling exponents x and C given by 

a..3) 

exactly, for all spatial dimensions d in the range 2 < d < A. For (i>4,x = l — | and C = 1) 
as in the isotropic case, while in c? = 2, where the model becomes identical to the isotropic 
model (the easy plane of motion being the entire space in that case), we again recover the 
isotropic results, C = f ? and x = — |. For the physical case d = 3, we have 
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X{d = 3) = -^ ; ad = 3) = ^ (1.25) 

The Fourier-transformed, equal time density (p — p) correlation function Cp{q) also obeys a 
scaling law 

CM = + ir'ft f Y^%,1 (1.26) 

where Ya{9xy) is a finite, non-zero, 0(1) function of the angle = ta.n^^ {q^ / qy) , and the 
scaling function follows 



constant, 



''i^^^^f (^^^^f) .e,{ql + .\qu?)»i^.q.f 

where ^Q is a "microscopic" length (of the order the interbird distance) and v a dimensionless 
nonuniversal constant of order unity. 

Finally, the hydrodynamic mode structure of this anisotropic flock consists of a pair of 
propagating longitudinal sound modes, with dispersion relation given, in the co-ordinate 
system of figure (1^) , by 



a; = c± (%, 0„-.) q - iqlfA ( , (1.28) 

where /a is a universal scaling function, c± {6^, (p^) is given by equation ( |6.22| ) of Chapter 6 
and the dynamical exponent 

z = 2C = = - (1.29) 

where the last equality holds in = 3. Note that this value of z again reduces to that of the 
isotropic model m d = 2, and d = A. 

The spatio-temporally Fourier transformed density-density correlation function 
(|p(g, C(j)p) has the same structure as that illustrated for the isotropic problem in figure 
(4), with the modification that is replaced by q^, and the scaling function is replaced 
by /a. The detailed expression for (|p(g, u;)p) is given by equation (|6.31|) of Section VI. 
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The remainder of this paper is organized (and we use the term loosely) as follows: In 
Section II, we formulate the isotropic model. In Section III, we specialize this model to 
the "broken symmetry" state, in which the flock is moving with a non-zero mean speed 
(v). In Section IV, we linearize the broken symmetry state model, and calculate the cor- 
relation functions and scaling laws in this linear approximation. In Section V, we study 
the anharmonic corrections in the broken symmetry state, show that they diverge in spatial 
dimensions d < 4, derive the new scaling laws that result in that case, calculate the exact 
exponents in d — 2 and discuss the difficulties that prevent us from obtaining these expo- 
nents for 2 < d < 4. In Section VI, we repeat all of the above for the anisotropic model. In 
Section VII, we describe in some detail how our predictions might be tested experimentally, 
both by observations of real flocks of living organisms, and in simulations. And finally, in 
Section VIII, we discuss some of the open questions remaining in this problem, and suggest 
some possible directions for future research. 



II. THE ISOTROPIC MODEL 

In this section, we formulate our model for isotropic flocks. As discussed in the intro- 
duction, the system we wish to model is any collection of a large number N of organisms 
(hereafter referred to as "birds") in a (/-dimensional space, with each organism seeking to 
move in the same direction as its immediate neighbors. 

We further assume that each organism has no "compass;" i. e., no intrinsically preferred 
direction in which it wishes to move. Rather, it is equally happy to move in any direction 
picked by its neighbors. However, the navigation of each organism is not perfect; it makes 
some errors in attempting to follow its neighbors. We consider the case in which these errors 
have zero mean; e. g., in two dimensions, a given bird is no more likely to err to the right 
than to the left of the direction picked by its neighbors. We also assume that these errors 
have no long temporal correlations; e. g., a bird that has erred to the right at time t is 
equally likely to err either left or right at a time t' much later than t. 
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Although the continuum model we propose here will describe the long distance behavior 
of an?/ flock satisfying the symmetry conditions we shall specify in a moment, it is instructive 
to first consider an explicit example: the automaton studied by Vicsek et al ||1|. In this 
discrete time model, a number of boids labeled by i in a two-dimensional plane with positions 
{rj(t)} at integer time t, each chooses the direction it will move on the next time step (taken 
to be of duration At = 1) by averaging the directions of motion of all of those birds within 
a circle of radius Rq (in the most convenient units of length i?o = 1) on the previous time 
step (updating is simultaneous). The distance Rq is assumed to be -C L, the size of the 
fiock. The direction the bird actually moves on the next time step differs from the above 
described direction by a random angle rii{t), with zero mean and standard deviation A. The 
distribution of T]i{t) is identical for all birds, time independent, and uncorrelated between 
different birds and different time steps. Each bird then, on the next time step, moves in the 
direction so chosen a distance voAt, where the speed vq is the same for all birds. 

To summarize, the rule for bird motion is 

e,it + l) = {e,it)) + v^{t) (2.1) 

fi{t + l) = fi{t) +vo{cose{t+ 1), sine (t + l)) (2.2) 

{r],{t)) = (2.3) 

{Vi{t)Vjit')) = A6,,6u' (2.4) 
where the average in eqn ( pTTP is over all birds j satisfying 

\rjit) - f,{t)\ < Ro (2.5) 

and 6i{t) is the angle of the direction of motion of the zth bird (relative to some fixed 
reference axis) on the time step that ends at t. The fiock evolves through the iteration of 
this rule. Note that the "neighbors" of a given bird may change on each time step, since 
birds do not, in general, move in exactly the same direction as their neighbors. 
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This model, though simple to simulate, is quite difficult to treat analytically. Our goal 
in our previous work |^ and this paper is to capture the essential physics of this model in 
a continuum, "hydrodynamic" description of the flock. Clearly, some short-ranged details 
must be lost in such a description. However, as in hydrodynamic descriptions of equilibrium 
systems Q, as well as many recent treatments |^ of non-equilibrium systems, our hope is 
that our continuum approach can correctly reproduce the long-distance, long-time properties 
of the class of systems we wish to study. This hope is justified by the notion of universality: 
all "microscopic models" (in our case, different specifications for the exact laws of motion 
for an individual bird) that have the same symmetries and conservation laws should have 
the same long distance behavior. This belief can be justified by our renormalization group 
treatment of the continuum model. 

So, given this lengthy preamble, what are the symmetries and conservation laws of flocks? 

The only symmetry of the model is rotation invariance: since the "birds" lack a compass, 
all direction of space are equivalent to other directions. Thus, the "hydrodynamic" equation 
of motion we write down cannot have built into it any special direction picked "a priori"; 
all directions must be spontaneously picked out by the motion and spatial structure of the 
flock. As we shall see, this symmetry severely restricts the allowed terms in the equation of 
motion. 

Note that the model does not have Galilean invariance: changing the velocities of all the 
birds by some constant boost {4 does not leave the model invariant. Indeed, such a boost 
is impossible in a model that strictly obeys Vicsek's rules, since the speeds of all the birds 
will not remain equal to fo after the boost. One could image relaxing this constraint on the 
speed, and allowing birds to occasionally speed up or slow down, while tending an average 
to move at speed Vq. Then the boost just described would be possible, but clearly would 
change the subsequent evolution of the flock. 

Another way to say this is that birds move through a resistive medium, which provides a 
special Galilean reference frame, in which the dynamics are particularly simple, and different 
from those in other reference frames. Since real organisms in flocks always move through 
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such a medium (birds through the air, fish through the sea, wildebeest through the arid dust 
of the Serengeti), this is a very reahstic feature of the model. 

As we shall see shortly, this lack of Galilean invariance allows terms in the hydrodynamic 
equations of birds that are not present in, e. g., the Navier-Stokes equations for a simple 
fluid, which must be Galilean invariant, due to the absence of a luminiferous ether. 

The sole conservation law for flocks is conservation of birds: we do not allow birds to be 
born or die "on the wing" . 

In contrast to the Navier-Stokes equation, there is no conservation of momentum. This 
is, ultimately, a consequence of the absence of Galilean invariance. 

Having established the symmetries and conservation laws constraining our model, we 
need now to identify the hydrodynamic variables. They are: the coarse grained bird velocity 
fleld v{r, t), and the coarse grained bird density p(r, t). The fleld v{r, t), which is deflned for 
all r, is a suitable weighted average of the velocities of the individual birds in some volume 
centered on r. This volume is big enough to contain enough birds to make the average 
well-behaved, but should have a spatial linear extent of no more than a few "microscopic" 
lengths (i. e., the interbird distance, or by a few times the interaction range Rq). By suitable 
weighting, we seek to make v{f, t) fairly smoothly varying in space. 

The density p{r,t) is similarly deflned, being just the number of particles in a coarse 
graining volume, divided by that volume. 

The exact prescription for the coarse graining should be unimportant, so long as p{r,t) 
is normalized so as to obey the "sum rule" that its integral over any macroscopic volume (i. 
e., any volume compared with the aforementioned microscopic lengths) be the total number 
of birds in that volume. Indeed, the coarse graining description just outlined is the way 
that one imagines, in principle, going over from a description of a simple fluid in terms of 
equations of motion for the individual constituent molecules to the continuum description 
of the Navier-Stokes equation. 

We will also follow the historical precedent of the Navier-Stokes equation by deriving 
our continuum, long wavelength description of the flock not by explicitly coarse graining the 
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microscopic dynamics (a very difficult procedure in practice), but, rather, by writing down 
the most general continuum equations of motion for v and p consistent with the symmetries 
and conservation laws of the problem. This approach allows us to bury our ignorance in a 
few phenomenological parameters, (e. g., the viscosity in the Navier-Stokes equation) whose 
numerical values will depend on the detailed microscopic rules of individual bird motion. 
What terms can be present in the EOM's, however, should depend only on symmetries and 
conservation laws, and not on the microscopic rules. 

To reduce the complexity of our equations of motion still further, we will perform a 
spatial-temporal gradient expansion, and keep only the lowest order terms in gradients 
and time derivatives of v and p. This is motivated and justified by our desire to consider 
only the long distance, long time properties of the flock. Higher order terms in the gradient 
expansion are "irrelevant" : they can lead to finite "renormalization" of the phenomenological 
parameters of the long wavelength theory, but cannot change the type of scaling of the 
allowed terms. 

With this lengthy preamble in mind, we now write down the equations of motion: 

dtv + Xi{v ■ V)v + X2{V ■ v)v + AsVd^/n = 

av-P\v\'^v - VP + £'bV(V ■ v) 

+DtV^v + D2{v-Vfv + f (2.6) 

oo 

P = P(p) = J2an{p-por (2.7) 

n=l 

1^ + V ■ (vp) = (2.8) 

where (3, Db, -D2 and Dj- are all positive, and a < in the disordered phase and a > in 
the ordered state (in mean field theory). The origin of the various terms is as follows: the 
A terms on the left hand side of eq. ( |2.6| ) are the analogs of the usual convective derivative 
of the coarse-grained velocity field v in the Navier-Stokes equation. Here the absence of 
Galilean invariance allows all three combinations of one spatial gradient and two velocities 
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that transform like vectors; if Galilean invariance did hold, it would force A2 = A3 = and 



Ai = 1. However, Galilean invariance does not hold, and so all three coefficients are non-zero 
phenomenological parameters whose non-universal values are determined by the microscopic 



rules. The a and (3 terms simply make the local v have a non-zero magnitude (= ^Ja/(3) in 
the ordered phase, where a > 0. -Dl,i,2 are the diffusion constants (or viscosities) reflecting 
the tendency of a localized fluctuation in the velocities to spread out because of the coupling 
between neighboring "birds" . The / term is a random driving force representing the noise. 
We assume it is Gaussian with white noise correlations: 

< Mr,t)f,ir\t') >= A6,,6\f-r')5it-t') (2.9) 

where A is a constant, and i , j denote Cartesian components. Finally, P is the pressure, 
which tends to maintain the local number density p(r) at its mean value po, and 6p = p — po- 



The flnal equation (|2.8| ) is just conservation of bird number (we don't allow our birds to 
reproduce or die "on the wing"). 

Symmetry allows any of the phenomenological coefficients Aj, a, an, P, Di in equations 
]6D and (p.7|) to be functions of the squared magnitude \vf of the velocity, and of the 



density p as well. 



III. THE BROKEN SYMMETRY STATE 

We are mainly interested in the symmetry broken phase; specifically in whether fluctu- 
ations around the symmetry broken ground state destroy it (as in the analogous phase of 
the 2D XY model). For a > 0, we can write the velocity field as v = foXy + 6v, where 
foa;|| =< V > is the spontaneous average value of v in the ordered phase. We will chose 
Vq = (which should be thought of as an implicit condition on vq, since a and /3 can, in 
general, depend on \vf); with this choice, the equation of motion for the fiuctuation 6v\\ of 
f II is 

dtSv\\ = —(Jid\\6p — 2a6v\\ + irrelevant terms. (3-1) 
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Note now that if we are interested in "hydrodynamic" modes, by which we mean modes 
for which frequency u; ^ as wave vector g ^ 0, we can, in the hydrodynamic {uj, q ^ 0) 
hmit, neglect dt6v\\ relative to a6v\\ in ( |3.1| ). The resultant equation can trivially be solved 

for 6v\\: 

Svi\ = -Dpd\\5p (3.2) 



where we've defined another diffusion constant Dp = Inserting ( |3.2| ) in the equations of 
motion for v± and 6p, we obtain, neglecting "irrelevant" terms: 



dtv± + 7(9||{7j_ + Ai (v± ■ V_l) v± + A2 (V± ■ v±^ 



= -V±P + DbVi, [V± ■ v±) + DTViv± + D\\dlvi_ + f± (3.3) 

^ + PoVx ■ ^± + Vx ■ {v^6p) + vod\\6p = Dpdldp (3.4) 
where Dp, Db, Dt and D\\ = Dx + D2VQ are the diffusion constants, and we've defined 

7 = Aifo. (3.5) 



The pressure P continues to be given, as it always will, by equation ( p.7|) . 

/^From this point forward, we will treat the phenomenological parameters Aj , 7, and Di 



appearing in equations (p.3|) and ( |3.4| ) as constants, since they depend, in our original model 



]6[),only on the scalar quantities \vf and p(r), whose fluctuations in the broken symmetry 
state away from their mean values Vq and po ^^re small. Furthermore, these fluctuations lead 
only to "irrelevant" terms in the equations of motion. 

It should be emphasized here that, once non-linear fluctuation effects are included, the 
Vo in equation ( |3.4| ) will not be given by the "mean" velocity of the birds, in the sense of 

{v)^^p^ , (3.6) 

where N is the number of birds. This is because, in our continuum language. 



/ p(r,t)v(r,t)d'^r) 

{J p{r,t)d'^r} (p) 
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while f in equation (p.5|) is 

vo = \{vif,t))\ 

Once p fluctuates, so that p =< p > +Sp, the "mean" velocity of the birds 



(3.8) 



(v) 



(pv) 




(p) 





(3.9) 



(P) (P) 

which only = Vq = \{v)\ if the correlation function {6pv) = 0, which it will not, in general. 
For instance, one could easily imagine that denser regions of the flock might move faster; in 
which case {6pv) would be positive along (v) . Thus, (v) measured in a simulation by simply 
averaging the speed of all birds, as in equation ( |3.6| ), will not be equal to vq in equation 
( 3.5|) . Indeed, we can think of no simple way to measure Vq, and so chose instead to think 
of it as an additional phenomenological parameter in the broken symmetry state equations 
of motion ( p.3|) . It should, in simulations and experiments, be determined by fitting the 
correlation functions we will calculate in the next section. One should not expect it to be 
given by {v) as defined in equation ( p.7|) . 

Similar considerations apply to 7: it should also be thought of as an independent, phe- 
nomenological parameter, not necessarily determined by the mean velocity and non-linear 
parameter Ai through 



IV. LINEARIZED THEORY OF THE BROKEN SYMMETRY STATE 



As a first step towards understanding the implications of these equations of motion, we 
linearize them in v± and 6p = p — p^. Doing this, and Fourier transforming in space and 
time, we obtain the linear equations 



(uj - 7g||) + Ft (q)] Vt (g, = h (g, t^) 



i\uj- 7^11 j + Tl {q)\ Vl + i(Tiq±6p = fi {q, uj) 
i {u - VQq\^ + Fp(g) 5p + ipoq±VL = 



(4.1) 
(4.2) 
(4.3) 
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where 



vl {q, uj) 



q_L ■ v± {q, uj) 



Q± 



(4.4) 



and 



Vt {q, uj) = v± {q, uj) - 



q±VL 
<1± 



(4.5) 



are the longitudinal and transverse (to gl) pieces of the velocity, Jt {q,oj) and {q,uj) are 
the analogous pieces of the Fourier transformed random force f{q,u!), and we've defined 
wavevector dependent transverse, longitudinal, and p dampings I'l,t,p' 



Tl {q) = DLqi + D\\ql 



(4.6) 



Tt (q) = Drql + D^qf^ , 



(4.7) 



Tp (q) = Dpql , 



(4.8) 



where we've defined Dl = Dt + Db, q± = \q±\- 

Note that in d — 2, the transverse velocity vt does not exist: no vector can be per- 
pendicular to both the x\\ axis and q± in two dimensions. This leads to many important 
simplifications in d = 2, as we will see later; these simplifications make it (barely) possible 
to get exact exponents in d — 2 for the full, non-linear problem. 

The normal modes of these equations are d — 2 purely diffusive transverse modes associ- 
ated with Vt, all of which have the same eigenfrequency 



(4.9) 



and a pair of damped, propagating sound modes with complex (in both senses of the word) 
eigenfrequencies 



a;± = c± (%) q - ITl 



v± {9g 



.2c2 (9^) 

= c± {e^)q-t{DLql + D^ql) 



2c2 (9^) 



v± 



2c2 (9^) 



iDpqj 



2c2{9^) 



(4.10) 
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where 9^ is the angle between q and the direction of flock motion (i. e., the x\\ axis), 

c± {e^) = cos (%) ± C2 {Oi) (4.11) 

v± (%) = ±^-^ COS (%) + C2 (^,-) (4. 12) 



C2(%) = ^i(7-t^o)'cos2(%) + c2sin2(e,^) , (4.13) 

and Co = ^/o\fH). A polar plot of this highly anisotropic sound speed is given in flgure 2. We 
remind the reader that here and hereafter, we only keep the leading order terms in the long 
wave length limit, i. e., for small q\\ and q^. 

The linear equations ([4. 1|) - (|4.3| ) are easily solved for the flelds 5p, vt and fL in terms 
of the random forces: 

Vt (g, Lo) = Gtt (g, tu) /t (g, tu) (4.14) 

vl (g, u;) = Gll (g, ^) /l (g, uj) + G^^ (g, u;) fp (g, u;) (4.15) 

5p (g, uj) = GpL (g, uj) fi (g, ^) + Gpp (g, ^) /p (g, uj) (4.16) 
where the propagators are 

Gtt = —J \ (4.17) 

-t [u - 7g|| j + Ft (gj 

G.. = , ^(.-^,||)-r (a ^^^^^^ 



(a; - c+ (%) g) (^ - c_ 


(%)g)+za; (r^ {q) + Tp (g)) 
icrig_L 




{voTl {q)+irp (g)) 


(t^ - c+ (%) g) {uj -C-{ 


;^,-.)g)+zcu (r^ {q} + Tp (q)) 


-m 


(voTl iq)+lTp (q)) 


{u - c+ (%) g) (tu - c_ ( 






{voTl {q)+irp {q}) 



= .,;r(-,^ (4-20) 
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(uj - 7g||) - Tl (q) 

{uj - c+ (6^) q) (cj - c_ {e^) q) + iuj {T l {q)+Tp{q))- iq\\ {vqT l (q) + iTp (q)) 



In writing the definitions of tlie propagators ( [4.14|) -( p:.16|) , we liave introduced a fictitious 



force fp in tlie p equation of motion ([4.3|) . Of course, tliis force is, in fact, zero; but the 
propagators Gpp and Glp nonetheless prove useful in the perturbative treatment of the 
non-linear corrections to this linear theory, so we have included fp here. 

Given the expressions ( [4.14| )-( P:.21| ) for the velocity and density in terms of the random 
force /, and the autocorrelation ( |2.9| ) of that random force, it is straightforward to calculate 
the correlations of the densities and velocities. We find: 

dj (g, uj) = (^v^ {-q, -u) vf (g, u)) 

= Gtt {q, ^) Gtt (-g, -w) (M ^) /t, {~q, 



+ Gll (q, t^) Gll {-q, -t^) ^^-Y- ifL (<?, ^) fh {-q, -uj)) 

Gtt (g, uj) (g) + Gll (g, uj) Lf^ (g) (4.22) 



where 



Lf, (g) = ^ (4.23) 

gi 

i^^ (g-) = 5j - Li^ (g-) (4.24) 

are longitudinal and transverse projection operators that project any vector perpendicular 
to both the flock motion and q±, 

Gtt (q, uj) = ^ (4.25) 

(c^-7g|i) +m<i) 

and 

Cll {q,uj) = -— ^ ^ (4.26) 

{uj - c+ (e^) qf {uj - c [e^) qf + (u {T l (q) + (q)) - gy {voTl (q) + iTp (q)) 

The transverse and longitudinal correlation functions Eqn. ( [4.25] ) and ( [4.2(j| ) are plotted 
as functions of uj for fixed q in figure 6. Note that they have weight in entirely different 
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regions of frequency: Ctt is peaked at = 7g||, while Cll has two peaks, dX u = c± (9^) q. 
Since all three peaks have widths of order g^, there is little overlap between the transverse 
and the longitudinal peaks as |g| ^0. 

The density-density correlation function 



C, 



pp 



{u - c+ (%) qf {u - c_ (%) qf + (u {Tl [q) + Tp (q)) - q\\ [voTl (q) + jTp [q 
looks almost identical to Cll, especially when one notes that near the frequencies u 



(4.27) 



(9,) q where both peak, the numerator of (= (c, (9,) q - .o.„)^) , differs from that 

of ( |4.27| ) only by a \q\ independent factor of (c± (%) q — fo?||) /poO'iQ'i- 

Given these Fourier transformed correlation functions, it is straightforward, and instruc- 
tive, to Fourier transform back to real time. In particular, it is simple to calculate the 
spatially Fourier transformed equal time velocity correlation function: 

duj „ ,^ . ^ t , ^ du 



^Ctt {q, UJ) + Lj-^ (q) / ^Cll {q, u;) 



A 
2" 



oc 



(4.28) 



TTiq) "^'^'TLiq), 

where the second integral over frequency has been evaluated in the limit of |g| ^ 0, so 
that c{9^) q ^ Tl oc q^, and the factor (p (g) depends only on the direction q of q, not its 
magnitude, and is given by the sadly complicated expression 



iq) = 



C2 (9^) q 



c+ {9q) q - VQq\\ 



+ 



c+ (%) q - Voq\\ + (c+ {9^) q - 7g||) ^ 
c- (%)?-^og||)^ 



c_ {9^) q - VQq\\ + (c_ {9^) q - 7g||) f 



1 



F (g, K, 7) 
+ 



A\ (g, 7) 



(g-',/€,7)-A_ (g-',«:,7)l£i| 
^-('?,«^,7) 



A_ (g-',K,7)-A+ (g,«:,7)I^ 



(4.29) 



where we've defined 



F(g;K,7) 



7-^0 



1 



(4.30) 
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{q- n, 7) = ±F (g; n, 7) + (\^) , (4-31) 



and 



(4.32) 



The second equality in (|4.29|) is obtained from the first simply by cancelling common factors 



of cTipoQ'l out of the numerator and denominator of various terms. 

Note, and this will prove to be crucial later, that (p {q) depends only on q, diffusion 
constants, and the dimensionless ratios k and 7/^0- This last fact is essential for our renor- 
malization group scaling analysis, as we will show later. 

The ^ divergence of ( [4.28| ) as |g| — >• reflects the enormous long wavelength fluctuations 
in this system. 

These fluctuations predicted by the linearized theory are strong enough to destroy long 
ranged order in d < 2. To see this, calculate the mean squared fluctuations in v\_ {f,t) at a 
given point r, and time t. This is simply the integral of the trace of Eqn ( ^.28| ) over all q: 

{\v±ir,t)\'^) = j -^^{vi{q,t)vi{-q,-t)) 

(4.33) 



A 



d'^q 
("2^ 

d'^q I id-2) (j){q) 



{27iy \DTqi + D\\qf^ DLqi + D\\qf^ 
The last integral clearly diverges in the infrared (|g| 0) for d < 2. The divergence in the 
ultraviolet (|^ ^ 00) for c? > 2 is not a concern, since we don't expect our theory to apply 
for |g| larger than the inverse of a microscopic length (such as the interaction range io). 
Presumably, at larger wavenumbers, the correlation function falls off fast enough that the 
wavevector integral in Eqn (|4.33|) converges in the ultraviolet. 

Indeed, we will in subsequent calculations mimic the effect of this putative more rapid 
decay of correlations as |g| — 00 with a sharp ultraviolet cutoff. We will restrict integrals 
over wavevectors to hypercylindrical shell with long {very long!) axis along the direction of 
flock motion xf 

\q±\ < A, -00 < gii < 00 (4.34) 
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with the ultraviolet cutoff A of order the inverse of a microscopic length (e. g., io)- 

Obviously, this is quite an arbitrary choice of ultraviolet cutoff, and any result that 
depends on the precise form of this cutoff will not be accurately calculated by this prescrip- 
tion. However, universal, long wavelength properties of the flock should be unaffected by 
the precise choice of cutoff, and it is on those properties that we will focus our attention. 

The infra-red divergence in Eqn ( |4.33|) for d < 2 cannot be dismissed so easily, since 
our hydrodynamic theory should get better as |g| 0. Indeed, in the absence of non- 
linear effects, this divergence is real, and signifies the destruction of long ranged order in the 
linearized model by fluctuations, even for arbitrarily small noise A, in spatial dimensions 
d < 2, and in particular in d = 2, where the integral in Eqn ( |4.28| ) diverges logarithmically 
in the infra-red. This is so since, if (^l"?!!^^ is arbitrarily large even for arbitrarily small A, 
our original assumption that v can be written as a mean value (v) plus a small fluctuation 
v± is clearly mistaken; indeed, the divergence of v± suggests that the velocity can swing 
through all possible directions, implying that (v) = for d < 2. 

In d = 2, this result is very reminiscent of the familiar Mermin-Wagner-Hohenberg 
(MWH) theorem which states that in equilibrium, a spontaneously broken continuous 
symmetry is impossible in d = 2 spatial dimensions, precisely because of the type of loga- 
rithmic divergence of fluctuations that we have just found here. 

In the next section, we will show that this prediction is invalidated by non-linear effects, 
and, in fact, much of the scaling of correlation functions and propagators is changed from 
that predicted by the linearized theory in spatial dimensions d < 4. 



V. NON-LINEAR EFFECTS AND BREAKDOWN OF LINEAR 
HYDRODYNAMICS IN THE BROKEN SYMMETRY STATE 
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A. Scaling analysis 



In this section we analyze the effect of the non-hnearities in equations (|3.3| ) and 
on the long length and time behavior of the system, for spatial dimensions d < 4. We will 
rescale lengths, time, and the fields v± and 6p according to 

x± — *• bx± 

x\\ — > b'^x\\ 

t^bH 
v± — > b'^v± 

5p^h^'5p (5.1) 

choosing the scaling exponents to keep the diffusion constants Db,t,p,\\i and the strength A 
of the noise fixed. The reason for choosing to keep these particular parameters fixed rather 
than, e.g., di, is that these parameters completely determine the size of the equal time 
fluctuations in the linearized theory, as can be seen from Eqn. ( [4.33|) . Under the rescalings 



the diffusion constants rescale according to Db,t V ^Db,t and Dp^\\ — > b^ ^^D 



hence, to keep them fixed, we must choose z = 2 and Q = 1. The rescaling of the random 
force / can then be obtained from the form of the f — f correlations Eqn. ( ^2.91 ) and is, for 
this choice of z and (: 

f^b-'-'/'f (5.2) 



To maintain the balance between / and the linear terms in v± in eq. ( |3.3|) , we must rescale 
the velocity field according to: 

vi_ b^v± (5.3) 

with 

X = l-d/2 (5.4) 
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which is the roughness exponent for the hnearized modeL That is, we expect v± fluctuations 
on length scale L to scale like L^. Therefore, the linearized hydrodynamic equations, ne- 
glecting the nonlinear convective term and the non-linearities in the pressure, imply that v± 
fluctuations grow without bound (like L^) as L ^ oo for (i < 2, where the above expression 
for X becomes positive. Thus, this linearized theory predicts the loss of long range order in 
d < 2, as we saw in Section IV by explicitly evaluating the real space fluctuations. 

Making the rescalings as described in Eqns. the equation of motion (|3.3| ) becomes: 



dtv± + 6^"7(9||til + [Ai {vi_ ■ V±) vi_ + Aa (v± ■ v±j vi_ 

= ( E b-'-aniSpr] + DbV^Vi. ■ VI.) + DtVIv± + D^idps, + h (5-5) 



\n=l 



with 



7a = X + 1 = 2 - d/2 (5.6) 
7, = ^ - C = 1 (5.7) 



and 



7n = z-x + nx-l=n+ (l-n)^ . (5.8) 

The scaling exponent Xp ^oi 6p is given by Xp = since the density fluctuations 6p are 
comparable in magnitude to the v± fluctuations. To see this, note that the eigenmode of 
the linearized equations of motion that involves 6p is a sound mode, with dispersion relation 
u = c± (6^) q. Inserting this into the Fourier transform of the continuity equation (|3.4]) , we 
see that Sp ~ '^^^^ . The magnitude of q± drops out of the right hand side of this expression; 
hence Sp scales like at long distances. Therefore, we will choose Xp = X = 1 ~ f ■ 

The first two of these scaling exponents for the non-linearities to become positive as 
the spatial dimension d is decreased are 7a and 72, which both become positive for d < 
4, indicating that the Ai (jJ± ■ v±, A2 (v_l ■ v±j v± and (J2V ±{6p'^) non-linearities are all 
relevant perturbations for d < A. So, for d < A, the linearized hydrodynamics will break 
down. 
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What can we say about the behavior of Eq. (|3.3|) and ( |3.4]) for d < 4, when the hnearized 
hydrodynamics no longer holds? The standard approach for such problems is the dynamical 
renormalization group. In most cases, this approach is only practical near the upper critical 
dimension (in our case dc = 4), and yields the anomalous exponents in an expansion in 
e = dc — d. This approach will obviously not be of much use in our problem in d = 2, where 
the ostensibly small parameter in this expansion e = 2. We will nonetheless undertake this 
approach in subsection B, and show that, for unfortunate technical reasons, we learn little 
even near d = 4. Fortunately, as we show in subsection C, because of the various symmetries 
in Eq. ( |3.3|) , we can obtain the exact scaling exponents in = 2. 

B. Renormalization Group Analysis, d < 4 

In this subsection, we analyze the effect of the relevant non-linearities Ai, A2, and (T2 on 
the broken symmetry state in spatial dimensions d < 4. 

Our tool is the dynamical renormalization group (for details, see, e.g., the excellent 
description in Forster, Nelson and Stephen 0). We will summarize the essential features of 
this procedure here; readers interested in details are referred to P]. 

We proceed through the iteration of the following 3 steps: 

1. We separate the fields v± and 6p, and the random forces / into short and long wave- 
length components, according to: 

v^> ir, ^) = X 1^ i^^ -) e'''-'-^'' (5-9) 

v±< (r, t)= v± ^) e^^---*^ (5.10) 

J< (27r) 

where denotes a wavevector integral restricted to a hypercylindrical shell b^^A < 
|(fj_| < A, where A is an ultraviolet cutoff, and likewise denotes an integral over the 
interior of this shell: < b^^A. 6p and f± are likewise separated. 
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2. Average the EOM over the short wavelength fields v±y, Spy, and /> to get new, 
eff'ective EOM for the long- wavelength fields and 5p<, with "intermediate" renor- 
malized parameters D^^, etc. This average is performed perturbatively in the non- 
linearities in the EOM. The perturbation theory can be represented graphically; the 
interested reader is referred to the previously mentioned [^] for further details on the 
mechanics of this. 

3. We now rescale the time, space, and the fields in the EOM according to ( |5.1|) in order 
to restore the original ultraviolet cutoff A of the problem. We will choose rescaling 
exponents z, ( and x to produce fixed points. 

Of course, the exponents are, in fact, completely arbitrary. We need not choose them to 
produce fixed points. However, it is very convenient to do so, since, as we will show in more 
detail later, the values of z, (, and x that do produce fixed points are exactly the values 
of the physical observable time, anisotropy, and roughness exponents that characterize the 
scaling properties of various correlation functions. 

Performing this RG procedure, we find the following recursion relations: 



^^ = {z-2 + G^^A9))Db,t (5.11) 

^ = (^-2C + G||,p(te}))^||,P (5.12) 

^ = {z + {n-l)x-l + G: {{g,})) (5.13) 

^ = (--l)po (5.14) 

^ = {x-l + z + {{g^})) Ai,v (5.15) 

^ = (^ - C - 2x + 1 - + Ga {{g^})) A (5.16) 
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^^(-0^.0 (5.17) 

§ = (^-07 (5.18) 

where we've taken b = 1 + di, di <^ 1, to obtain differential recursion relations, the G's 
represent graphical (i. e., perturbative) corrections, and the {giYs are a set of dimensionless 
coupling constants involving ratios of powers of the dynamical parameters. We have also 
dropped "irrelevant" terms in these recursion relations. The coupling constant Xp is the 
coefficient of the V_l ■ {v±6p) non-linearity in the p equation of motion. This coupling 
constant is equal to 1, and must, up to trivial rescaling corrections, remain equal to 1 
upon renormalization. This is a simple consequence of the fact that mass conservation is 
exact; that is, the equation of motion for p must remain the simple continuity equation 
dtp + V ■ (pv) = 0, except for the trivial changes introduced by rescaling. This implies that 
Gp {{gi}) = 0, exactly, for all {gi}. We will later use this fact to obtain exact values for the 
scahng exponents z and ( in d = 2. 

It is worth noting that 7 is treated as an independent variable here, only its bare value 
is related to the bare values of Ai and fo through eq. (p.5|). 

Although there is no symmetry argument forbidding renormalization of fo and 7, simple 
power counting shows that there are no relevant graphical corrections to them; this is why 



no graphical corrections appear in (|5.17|) , (|5.18|) . 



As mentioned earlier, the rescaling exponents z and ( are arbitrary. We chose them 
to produce fixed points only for computational convenience. However, there is no choice of 
rescaling exponents that will keep all the parameters fixed. For instance, to keep Db,t and 



D\ip fixed, we will have to choose z > 1. However, with this choice of z, the relations (|5.13|) 



and ( |5.14| ) show that ai and po flow to infinity. 

Which of the parameters, then, should we choose keep fixed? That is, which is most 
convenient to keep fixed? The answer to this question is provided by the renormalization 
group matching formalism. This approach enables one to use the renormalization group to 
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relate correlation functions in the original, unrenormalized model at long distances and large 
times to the same correlation functions in the renormalized system at shorter distances and 
times. The advantage of this approach is that long distance, large time correlation functions 
are hard to calculate in d < 4, since, as we showed from our earlier scaling arguments 



are not accurately calculable from the harmonic theory developed in section IV, since the 
non-linearities Ai,2, Ap, and (J2 have very large effects at long distances. By mapping these 
correlation functions onto those at short distances in the renormalized equations of motion, 
we circumvent this problem. Clearly, there is a caveat here: even at short distances, the 
correlation functions in the renormalized model can only be calculated accurately in the 
harmonic theory if the non-linear couplings in that renormalized model are not too big. 
This suggests that the convenient choice of the rescaling exponents z and ( is that which 
keeps the non-hnearities Ai,2, Ap, and a2 fixed. 

Let us illustrate these considerations explicitly for one very important correlation func- 
tion: the equal time spatially fourier transformed velocity-velocity autocorrelation function: 



This particular correlation function is important because it gives us our best measure of 
the size of the velocity fluctuations, and will ultimately determine whether or not these 
fluctuations destroy the long ranged orientational order of the flock (thereby driving its 
mean velocity to zero). 

Cij (q) is, of course, a function of the flock dynamical parameters Db,Ti -D||,p, A, etc., 
as well as of q. Furthermore, at small g, it is difficult to calculate in spatial dimension 
d < A due to the non-linear terms, for the reasons discussed above. So let's follow this 
renormalization group matching procedure to relate Cij {q; {B^}) where {B^} denotes the 
set of dynamical parameters D%rp, Dp Aq, etc. in the unrenormalized model, to the same 
correlation function in the renormalized model, a renormalization group time i later: 



(and can also verify from the renormalization group recursion relations ( |5.11| )-( [5.12| )), these 




(5.19) 




(5.20) 
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where the {Bi{i)} denote the renormahzed parameters. 

In the discussion that follows, we will first consider the case (^^^ <^ (^)^ ' 
conclusion of the discussion of this special case, we will briefly indicate how the general 
case can be treated to obtain the scaling laws quoted in the introduction. For the case 
( A^) ^ (^)''' ^^^^ choose i = {q±) = In (^-^^, where A is the ultraviolet cutoff, on the 
right hand side, and obtain 

/ A \ 2x+C+rf-i fa \ 
= (^-j C,, \ a, {B^L (q^))} j (5.21) 

Now, 2/ the original q± was small (^ A) and we have chosen the rescaling exponents xX^ 
and z so that the non-linearities Xi{i), A2(^), Ap(£) and o"2(£) on the right hand side of (|5.21|) 
flow, as £ — >• cx), to 0(1) fixed point values (A]; 2.p, '^2}, then Ai (£*), A2 (^*), Xp (^*) and (T2 (^*) 
can be replaced by those fixed point values, since will be large. Because those fixed point 
values are, by assumption, 0(1), then, up to 0(1) correction factors coming from these non- 



linearities, the right hand side of ( |5.21| ) can be evaluated in the harmonic approximation 
Eqns. ( |4.22| ), ( [4.25| ), ( [4.26|) . (The correction factors are only of 0(1) - i. e., not divergent - 
because the right hand side of ( p.21| ) is evaluated at large (f (|gl| = A), where the infrared 
divergences associated with the strong relevance of the non-linearities do not matter. It is 
precisely because of those infrared divergences that we could not evaluate the left hand side of 
( p.21| ) directly, but rather were forced to go through this seemingly circuitous RG matching 
formalism). Making that harmonic approximation on the right hand side, we obtain: 

\ (^) / d^a^ + d; 

I 



Cij 



A,^;{T{Q}\= -^Piiqi 



A, 

+ 



A,^;{B,iQ}\Li^iq^) (5.22) 



A 



where we have used the fact that we've chosen the scaling exponents to make A and all 
the diffusion coefficients {-Dj} flow to fixed points A*, {D*}. We wish to show that this 
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expression depends on q only through the scahng ratio u = 



. The first (transverse) 



term in ( 5.22 ) exphcitly has this property. The second (longitudinal term) would also, except 
for the (f) factor, to which we now turn. From Eqn. ( [4.29 ) for 0, we see that to calculate 
this factor we must calculate 
/ 



A 



\ 



( 



In 



'0 I** J 



(£* 



(5.23) 



A 



-AB^ (4)} 



A, 



(7 {L)-v^ (4)) 



tAB^ (4)} 

^^/s:(£. 



and 



B 



± 



A,-%;{i?,(4)}| =±F|A, 



\ (^) 



{5,(4)} 



, (7(t)-t;o(4)) 



(5.24) 



(5.25) 



all of which are clearly dependent only on the fixed point value of the ratio (which is just 
a number of 0(1) since 7 (£*) and vq (£*) have the same dependence on : exj) [(z — Q 
as can be seen from their recursion relations), and the combination 



(5.26) 



By combining the recursion relations ( |5.17| ), ( |5.13| ), and ( |5.14| ) into a recursion relation 
for K 



d , dlnvQ 1 (i , , ^ . 

^(/-'.) = ^-2^(/-a, + /npo) = l-C 



(5.27) 



we find 
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(5.28) 



which imphes that 



i-C / A ^ 



Using this in (|5.26|), we see that the combination: 



k{Q ^\ =ko^ . (5.30) 



Q± 



takes on precisely the value it would take on using the unrenormalized parameters and the 
unrescaled wavevector q. Hence, the same is true of F, A±, and B±. And, therefore, the 
same is true of (g). 

Thus, we can replace 0(A, ^ ; [Bi (£*)}) in (|5.22|) with (g; {5°}), its unrenormalized 



A 



value straight from the linearized theory. Doing so, and recalling that the unrenormalized 
(j) {q) was 0(1) for all directions q of q, we see from ( |5.22|) that the correlation function Cij 
is largest when 

For |g| < A, where our theory applies, ( [5. 311 ) implies that ^ q± (since ( < 1). In that 
limit, (f){q) ^ 1; using this in the expression ( |5.22D for Cij in the renormalized system, and 
using eq. ( p.22| ) in turn in our expression ( 5.21|) for Cij in the original model, we obtain, for 



gy ^ g±, the scaling law: 

I (^) J 



(5.32) 



Note that the range q\\ q± for which this scaling law holds includes those g's which 
dominate the fluctations; namely, those with ^ ~ (^)^" 

Integrating Cij (g) over all q gives the equal time, root-mean-squared real space fluctua- 
tion of v±: 

(K(?,*)0=/^a,® 
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A 



2x+d-l 

q± 



(2x+o!+C-l) 



dq\\ fi, 



(5.33) 



where the final proportionahty was obtained by scahng q± out of the q\\ integral via 



the change of variables q\\ = AQy 



and we've defined the g^-independent constant 



A = JdQiifij (Qii) ASx+'iK. The final integral in (|5.33|) clearly converges in the infra-red 
{\q±\ — > 0) limit if and only if x < 0. Furthermore, if x is > 0, and we impose an infra-red 
cutoff > Lj_'^ in ( |5.33| ), where is the lateral (i.e., _L direction) spatial extent of the 
system, we easily obtain 



(5.34) 



Indeed, the connected real-space, equal time, velocity autocorrelation function discussed 
in the introduction is given by 

Cc{r) = {v{r + R,t)-v{r,t)) - \{v{r,t))f 

[v± (^f + R,t) ■ v_i {f,t)] 
d'^q 



(27r)' 



(2vr) 

Making the changes of variable 

^1 = 



(5.35) 



Q- 



Q\\ 



Ri 



(5.36) 



we obtain the scaling law (|1.8|) for Cc (^Rj quoted in the introduction, with 



fv{u) 



■R±+Q 



•(2x+d+C-l) 



(5.37) 



This shows that the x we obtain from the renormalization group by the prescription 
we have chosen - namely, making the specific set of parameters Db,t,\\, -^i,2,p) and a2 
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flow to fixed points - is precisely the physical roughness exponent defined by the velocity 
fluctuations. 

To summarize: if we chose the rescaling exponents z, and ( so as to make the particular 
subset of the dynamical parameters Db,t,\\,p, A, Xi,2.p, and o"2 flow to non-zero fixed points, 
then those x and ( are the ones that appear in the scaling law ( |1.8| ). This directly and 
simply relates the RG to physically observable correlation functions, so this is the choice we 
will make. 

A scaling law similar to ( |5.32|) can be derived, by precisely the same type of arguments, 
for the equal-time density-density correlation function: 

CM - (^) YiM (5^38) 

where, in writing this relation, we have used the fact that the "roughness" exponent for p, 
Xp = Xi the "roughness" exponent for v^. 

The alert reader will have noticed that neither of the scaling laws ( |5.32| ) and ( |5.38| ) derived 
so far involve the time rescaling exponent z. This is unsurprising, since we've considered 
only equal time correlation functions up to now. 

To fully study the dynamics of the model, we need to consider correlations between 
different times, as well as positions. These different time correlation functions will involve 
z. 

It is easiest to work with the space and time Fourier transform Cij {q,uj), defined by 

{vi (g, uj) Vj (-g, u')) = 6 {u + u') Cij (g, u) (5.39) 

We will find, as we've asserted many times already in this paper, that Cij {q, u) does not 
have a simple scaling form, unlike the equal time correlations. Nonetheless we can derive an 
expression for it in terms of functions of q which do show simple scaling behavior; namely, 
effective wavevector dependent diffusion constants that diverge as |g| ,ti; 0. 

We begin this derivation by separating Cij into its transverse and longitude parts: 

a, (g, to) ^ LjjCL {q, CO) + P.^Ct {q, co) (5.40) 
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where Li 



1^ 1j 



— and = 6ij — L^- — are, respectively, the longitudinal and 

transverse projection operators defined in section IV. Both the transverse and longitudinal 
pieces Cl^t {q, ^) obey the same renormalization group transformation 



C,,T {q^,q\\,u;; {b^}) = e^--^^^^^'-'^'C,,T {/q^,e%,e^'u;; {B, (£)} 



(5.41) 



where we've been careful to take into account the rescaling of the delta function in ( p.39| ) in 
deriving the argument of the exponential in (|5.41|) . 

As for the equal-time correlation function, it is convenient here to chose the rescaling 
factor such that e^q± is right on the Brillouin zone boundary; i.e., = A. Making 

this choice, and taking A = 1, we obtain 



Cl,t {q±,q\\,uj; {Si}) = q^ 



g., ^, 4; {B. {Q} 



q± ?± 



(5.42) 



where, on the right hand side, we've defined 



= In 



\q±\ 



(5.43) 



For the moment, let's focus on the longitudinal piece Cl {q, uj). As we argued for the equal- 
time correlation function, here too we can evaluate the right hand side in the harmonic 
approximation Eqn. ( [4.26| ). This gives 



Ci(g-l,g||,^;{5°}) 



DEN 



(5.44) 



where 



DEN 



4 - ^+ (^*) 

ql 



4 - ^- iQ 
q± 



+ 



4 ) (r^ {Q + r, (Q) - vo {Q (tl {Q + ^^r, 
\qlj V ^o(4) 



(5.45) 



where 



uj± {Q = uj± ^g_L, 7 {Q , po (^*) , {Bi , 
with the sound frequencies uj± (g; 7, po, {Bi}) obtained in the harmonic theory: 



(5.46) 
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^± (g;7,Po, {^i}) = ^ ^ ^° g cos gg- ± 



\ 



(7 -"^0)9 COS gg- ^ 2-2/) 



(7 + ''^o) gii 



± 



'(7 - ^o)gii 



o"iPo91, 



(5.47) 



Si 



(5.48) 



(5.49) 



and D^, Dp and A* are the fixed point values of Dl, -Dy, and A, to wliicli tliose parameters 
will have flown for large, as it will be for small q±. 

The complication of this expression - that is, the fact that stops it from having a simple 
scaling form - is that the parameters 7 (^*), <Ji (i*), and po (^*) that appear implicitly in 
( ^.44| ) do not flow to field point values for our "canonical" choice of x, z, and (, as discussed 
earlier. Physically, this reflects the fact that the scaling of the sound speeds {u oc q) is 
different from that of their dampings (damping rate oc in harmonic theory; we will show 
damping rate is oc here, in a moment). 

To proceed, it is first useful to reorganize ( |5.44|) slightly; by multiplying numerator and 
denominator by q'^: 



CL{q±,qii,uj;{B^}) 



A^iuj - VoiQqiiqr')WI 



z-C\2jz-(;-2x+l-d) 



(5. 



Next, we solve the recursion relations for 7 (£*), ai (£=„) and po i^*)'- 



(4) = e^^-<^''Vo{i = 0) = VoiO)qi 



C.-Z 



(5.51) 



7(t) = e(^-^)^-7(^ = 0)=7(0)gi-^ 



(5.52) 



ai (4) = e(^-i)'^*ai(£ = 0)=ai(0)gi- 



1-2 



(5.53) 
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Po (t) = e(^-i)^*po(^ = 0) = po(0)g. 



l-z 



(5.54) 



where in the second equahty in each equation we've used = In ( • Using these results 
and the expressions ( p.46[ ) and ( p. 47] ) for uj±{i^), we see that z and C drop out of the 
combinations 



2 q\_ 



^ (7(0)-.o(0)) ^^C-A +.,(o)po(0)gf-) 



^ ^ , A7(0)-^^o(0))9|| 
(<?1, gy; 7(0), (0) , cTi (0) , Po (0) 



+ ai (0)po(0)gi 



(5.55) 



and therefore the positions of the peaks in the full correlation function (|5.5CI| ) are exactly 
those given by the harmonic theory using the bare parameters cri(O), fo(0), 7(0), and po(0), 
namely, w± (g; 7(0), fo(0), cri(O), po(0)). This is a direct consequence of the fact that there 
are no (relevant) graphical renormalizations of the parameters (7, cti, and po) that determine 
the sound speeds (see the recursion relations (|5.14|) - (|5.17|) ) and shows that the relevant non- 
linearities below d = 4 do not alter the positions of the peaks in the spatio-temporally 
Fourier-transformed velocity-velocity autocorrelations. 

The same cannot, however, be said for their widths. Indeed, using the above result for 
the sound speeds and Eqns. (|5.48| ) and ( p.49| ), we see that Cl {q,^^) can be rewritten: 



CLiq,uj) 



A,(a;-t;o(0)g||)2gr 



2jz-(;~2x+l-d) 



(5. 



(a; - c+ (9^) qf - c_ {6^) qf + [u: (rf (g) + (q)) - g,, {v^m^ (q) + 7(0)r« (g) 



where the sound speeds are given by the harmonic result equation ( [4.11j ), and the renormal- 
ized dampings 



rf(g1 



Si 



A. 



(5.57) 



r?(g) 



ql = qlfp 



(5.58) 
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obey simple scaling laws. 

The exact form of the scaling laws that we have obtained here (namely, e.g., fi 



D*^ y^j + ) , is not correct, because our choice of = ln{^) is only valid when q]_ 3> 
q\\. In the opposite limit qS^ q\\ , the fluctuations become negligible in the renormalized 
problem once -D||g^ becomes ^ DbA.^ in the renormalized problem, because at this point 
the linearized approximation to the correlation functions is smaller than its largest value at 
the Brillouin zone boundary. This means we can now stop the renormalization at i^, such 
that e'^^'qw = A x D*b/D^^ , which implies that = — ^ + 0(1), where the 0(1) factor 
is universal, because it depends only on the fixed point values of the diffusion constants. 
Performing the above calculations with this choice of i^, we now obtain 



r?(g1 













[4). 



4 X 0(1) . qlU (I 



(5.59) 



* / 1 



2-z-i 

.c \ — 



X 0(1). Note that the precise 



where we've now defined fLy^ 

form of this scaling function is different in this regime than that found earlier for q'^_i ^ q\\- 
Furthermore, its exact form is uncertain, due to our uncertainty in the 0(1) factor, which, 
as discussed earlier, cannot be determined without knowing the fixed point values D* of 
the diffusion constants. However, regardless of their values, we still get a scaling law with 
the same power of q± and the same scaling variable ^ as that found earlier in the opposite 
limit. 

In between these two limits we have to choose i^, to smoothly interpolate between the 
two limits. This choice will clearly depend on the ratio Naively, one could imagine 
simply choosing = ln{min{ — , ^t^)). A subtler choice would take into account the 0(1) 
perturbative corrections we've neglected, and would presumably lead to a smooth crossover 
of between the two limits. 

The moral of this discussion is three-fold: 

1) we always get scaling laws of the form ( |5.57| ) for the dampings, 

2) the renormalized damping functions and the noise strength are always of such a form 
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that they depend only on q\\ for q\\ ^ g^, and only on q± in the opposite hmit, and 

3) we can only calculate the scaling function if we know the diffusion constants at the 
fixed point. 

This last point will stop us from calculating the crossover functions ind = 2, even though, 
as we will see, we can calculate the exponents there. 

We see from ( |5.57] ) that the physical significance of the exponent z is that it gives the 
scaling of the peak widths (in uj) of Cl {q,uj) with q±, while the peak positions continue to 
obey the "z = 1" scaling u cc q. 

Similar, but actually far simpler, arguments show that the transverse correlation function 
Ct {q, uj) obeys 

{uj --fqiiY + T^{q) 

where Tt {q) obeys the scaling law 

Tt (q) = qlh (^^) , (5.61) 

and /a is a scaling function associated with A. For general values of q\_ and q\\, the same 
scaling function /a should also be present in all of the other correlation functions as well 
(wherever A appears), such as eq. ( [5.56| ) for the longitudinal correlation function. 



Likewise, the propogators of the full non-linear theory are given, in c? < 4, by the 
harmonic expressions (4.18)-(4.21), except that F^, Fp, and F^ in those expressions are 
replaced by the anharmonic scaling laws ( |5.57| ), ( |5.58| ) and ( p.61| ). 



To complete the specification of the scaling laws, we need the asymptotic behavior of the 
scaling functions /a,l,t,p(m)- From ( |5.57D , ( [5. 58] ), and the analogous result for ( |5.6lD , and 



requiring that the second point of our tripartite moral applies, we see that 

, constant, m — 
/aM«<I (5-62) 



u i , u ^ oo 



and 
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constant, m — > 

/l,t»cx<; ^ (5.63) 



which imphes that 



The simplest summary of the scahng of all correlation functions and propagators is: simply 
use the harmonic expressions for them, except that diffusion constants Dt,b,p should be 
replaced by wavevector dependent quantities that diverge as g — >• 0, according to the scaling 
law 

Dt,b,p (q) = ql-'fT,B,, I I , (5.65) 
the bare noise strength A should be replaced by 



A 

and the diffusion constant should be replaced by 

( 

V 



^11 {q) = iT^'h 

as can be seen by requiring that 



(5.67) 



D\\{q)ql = Dll^^^ ql (5.68) 

the right hand side being the form of the q\\ dependent term in ( ^.57) ). 

We hope the reader has not been too confused by the fact that we have restored the 
ultraviolent cutoff A ~ l/£o to the problem by going back to dimensionful units where 
A^l. 

This completes our discussion of how the renormalization group, and, in particular, the 
exponents Xi ^i^nd relate to physically observable correlation functions and propagators. 
Now, we turn to the problem of actually calculating those exponents. 
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C. Exponents in d = 2 



To do this, we must calculate the graphical corrections in equations ( 5.11 )-( 5.16|) . The 



procedure for this, as discussed in involves the harmonic correlation functions and prop- 
agators and vertices representing the non-linearities Xi,2,p and a2- Rather than actually 
calculating these corrections, we will show that, when A2 = 0, the structure of the theory is 
such that we can determine the exponents z and ( exactly. 

Consider first the Ai non-linearity. Separating v± into transverse and longitudinal com- 
ponents, 

= vt + vl (5.69) 

this can be written 

(v± ■ V^) v± = (vt ■ V_l) vt + (vl ■ V_l) Vl + (vt ■ V^) vl + (vl ■ V_l) vt (5.70) 

Now consider the graphs that can be constructed from vl — vt the cross-terms in this 
expression. These will always mix transverse and longitude propagators and correlation 
functions in the internal integrals over momentum and frequency. But, as noted earlier 
in our discussion of the harmonic theory, the peaks in the longitudinal propagators and 
correlation functions occur at different frequencies {uj = oj±{q)) than those in the transverse 
propagators and correlation function, which occur at = 0. Furthermore, the overlap 
between these peaks is negligible, since their widths (oc q\_ with ^ > 1) are much less than 
this offset in peak positions. This implies that the integral over wavevectors and frequencies 
of any graph that mixes transverse and longitudinal propagators and correlation functions 
will be much less (by powers of q) than any similar graph containing purely transverse or 
purely longitudinal propagators and correlation functions. Hence, the vl — vt cross terms 
in (|5.70|) are irrelevant compared to the pure vl and vt terms. 

Now let's consider those relevant pieces. The Fourier transform of the vt piece at wave 
vector q can be written in Fourier space: 
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FT ( [vt ■ V_l) VT,)^=iJ2^Tip} {q± - P±) vt, (q- P) 

p 

nf ^T, (p) Vt, (q-p) (5.71) 



p 



where we have used the fact that vt is transverse, so pj_ ■ vt (p) = 0. So this piece of the 
Ai vertex, which is a term in the equation for dtVi{q,t), is proportional to the external 
momentum q±. So is the purely longitudinal term, as can easily be seen in real space. Since 
vl is longitudinal, we can write 

vl = (5.72) 
for some scalar field 0. Now the second term in (|5.70|) can be rewritten in terms of 



{vL ■ Vx) vu = {d,<j)) {d,d,<l)) = {d,<j)) {d,dj<j)) = {d,<pd,<P) (5.73) 

which is clearly a total derivative, whose Fourier transform is proportional to gl. 

Hence, the two re/evant pieces of the Ai vertex are proportional to the external momentum 
q±. Clearly, the a2 term, being a total _L derivative, is also proportional to q± in Fourier 
space. Hence, when A2 = 0, all the remaining relevant vertices are proportional to q^. An 
immediate consequence of this is that A and D\\ acquire no graphical renormalization. For 
A, this can be seen by noting that any graph that renormalizes A (e.g., figure 7) must 
contain two external vertices each proportional to g^, and hence must be proportional to 
q]_. Therefore all renormalizations of A must be proportional to q]_, and hence negligible, 
as |g| 0, relative to the bare A. Likewise, any graph for the diffusion constants must be 
proportional to figure 8, which must be proportional to at least one power of q±. Since D\\ 
and Dp involve no powers of gl, they cannot he renormalized graphically. 

Thus, when A2 = 0, A, Dp and D\\ get no graphical renormalization. That is, Gy, Gp and 
Ga in (|5.12| ) and ( ^.16| ) are, exactly, = 0. Thus, the requirement that A, D\\, and Dp flow 



to fixed points i '^^l''' ^ ~ ^ ~ ^ leads to two independent exact scaling relations between 



di 



the three independent exponents x, z and C- Requiring = implies 

z = 2C (5.74) 
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while requiring ^ = leads to 



z = ( + 2x + d-l . (5.75) 



We emphasize that we have only shown that these relations ( |5.74|) and (|5.75| ) hold when 
A2 = 0. 

We can obtain a third independent exact scaling relation between these three exponents, 
and thereby determine them exactly, when A2 = 0, by considering the renormalization of 
the non-linearities Ai, Ap, and a2- 

We start deriving this third relation by noting that when A2 = and Ai = Ap = A, there 
can be no graphical renormalization of Ai. This is because for these parameter values the 
equations of motion ( |3.3| ) and ( p.4| ) have an exact symmetry which we call pseudo-Galilean 
invariance: namely, they remain unchanged under the "boost" transformation: 

f±^r± — Xvbt (5.76) 



v±ir, t) v±{r, t) + Vb (5.77) 

where the "boost" velocity Vb is an arbitrary constant vector in the _L plane. 

This symmetry must be preserved upon renormalization with the same value of A. Hence, 
there can be no graphical renormalization of A, when Ai = Ap. That is, = when 
Ai = Ap. Since = always, it is clear that, if Ai < Ap initially, it will always remain so 
upon renormalization. 

This implies that, for flocks that start with bare Ai less than the bare Ap (which should 
be a finite fraction of all possible flocks), only two types of fixed points are possible: 

I) Xl = X; = 0, or 

II) a; ^ 0. 

The first type of fixed point (I), however, is readily seen to be unstable to Ap, which must 
always be non-zero (and, in fact, = 1) initially. Hence, this fixed point is never reached, and 
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we must flow to a fixed point of type II. To see that fixed points of type I are unstable, note 
tliat for such a fixed point, the only remaining relevant non-linearity is (72 • But, by itself, a2 
can not renormalize any of the diffusion constants Db and Di. The reason for this is that 
any graph (e.g., figure 8) that renormalizes any diffusion constant must have an external 
velocity leg emerging from the right. However, using only the (T2 vertex, which only involves 
p, we can only make graphs with a p leg emerging from the right. Therefore, at a fixed point 
of type I, Gb,t = 0, exactly, in Thus, to find a fixed point for Db,t, we must choose 



z 



2. Combining this with the previous exact scaling relations ( ^.74|) and (|5.75|) , we find 



= 1 and X = 1 — |. But using these values (which are nothing but those that we found in 
the harmonic theory), we find that Ap is a relevant perturbation at any fixed point of type 
I: 

'-^^H)^' (5.78) 
for all spatial dimensions d < 4. Note that ( |5.78| ) is exact at any fixed point of type I, since 



Gp = 0, exactly. 

So fixed points of type I are unstable, and we will always flow to a flxed point of type II. 
Using the recursion relation (|5.15|) for Ap, and the fact that 0^ = always, we immediately 



obtain that A* can be 7^ if and only if a third exact scaling relation is satisfied, namely: 

X=l-z. (5.79) 



The three relations ( ^.74|) , ( ^.75|) , and ( p.79| ) that hold when A2 = can trivially be solved. 



to find the exact scahng exponents in all (i < 4 that describe fiocks with A2 = 0; 

d + 1 , , 

C = (5.80) 

2(d+l) 

z= ^ 7 ^ (5.81) 



and 



(5.82) 
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Note that these match continuously, at the upper critical dimension (1 = 4, onto their 
harmonic values ( = 1, z = 2, and x = l — f = — 1, as they should. 

If Xf"''^ > Xf"''^, then, besides the two cases that we discussed above, there may be a 
third type of fixed point: 

III) a; = 0, a* ^ 0. 

For this type of fixed point to be stable, the exponents x and z have to satisfy: 

X + z<l (5.83) 



The above inequality, together with eq. (|5.74| ) and eq. ( ^.75|) give eq. (|5.8CI|) , (|5.81| ) and 



( ^.82| ) with the "=" signs replaced by " < ". For simplicity, we will only discuss the cases 



with Af"'"'^ < Xp""^^, where the exponents are given by eq. (|5.80|) , (|5.81| ) and (|5.82|) . However, 



all the qualitative results will be valid for case III, e. g., the spontaneous symmetry broken 
phase will be more stable in case III if it is stable in case II because of the inequality. The 
simplest possible scenario is that there is only one stable fixed point, regardless of whether 
\Bare ^ \^are ^^^^ ^^^^ ^^ig,^ it is type II, and has the canonical exponents ( |5.80|) - ( p. 821) . 



We consider it highly probable that this is, in fact, the case. Even if it is not, ( |5.80| ) - ( |5.82| ) 
do hold for some flocks (those with A^"*"*^ > Af*^*"*^). 

Our derivation of these results ( |5.8CI| )-( |5782| ) depended only on the assumption that A2 = 
0. Note, however, that in c? = 2, any flock is equivalent to a flock with A2 = 0. This is 
because the Ai and A2 vertices become identical in c? = 2, where v± has only one component, 
which we'll take to be x. That is, in c? = 2, 

Ai (y^ ■ V^) vi_ = Xixv^d^v^ = ^Xid^ (vl^ x (5.84) 

A2 (V_L ■ v_l) vi_ = X2X (d^v^) = ^X2d^ (f^) X (5.85) 

so that the full v± non-linearity becomes | (Ai + A2) dx (f^) x, which is just what we'd get if 
we started with a (primed) model with Ag = and A'^^ = Ai + A2. This later model, since it 
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has A2 = 0, must have the "canonical" exponents ( |5.8CI|) - (|5.81|) hence, so must the (Ai,A2) 
model, which includes all possible d = 2 models. So all models in d = 2 must have the 
canonical exponents ( |5.80| )-( |5782D . 

Equivalently, we can derive this result by simply noting that, in d = 2, the full 
v_[_ — vertex becomes | (Ai + A2) ("W^.)) "which is a total x derivative even when A2 7^ 0. 
Furthermore, the d = 2 model now has the pseudo- Galilean invariance ( |5.76|) - (|5.77|) when 
Ai + A2 = Ap. These two properties {v_[_ — vertex total x derivative, and pseudo- Galilean 
invariance at a special point) are all that we used to derive the "canonical" exponent (|5.80|) - 
( ^.82| ); so those canonical exponents must hold in = 2. Setting = 2 in (|5.8CI| )- (|5.82|) , we 
obtain: 



c 



(5.86) 



(5.87) 



X 



1 

'5 



(5.88) 



Note, in particular, that x < 0. This implies, as discussed earlier, that the flock exhibits 
true long ranged order. 

Using the exponents (|5.86|) - (|5.88| ) in the general scaling relations, such as (5.31) and 
(5.33), we obtain all of the scaling results for correlation functions m d = 2 quoted in the 
introduction. Note also that for this set of exponents z — C, — 2x + — d = 0. Hence, 
from equation ( [5. 661 ), we see that the noise strength A is a constant, independent of q, 
which makes sense since A is unrenormalized graphically. So, in the d = 2 model, we can 
calculate all correlation functions from their harmonic expressions, except that we replace 
the diffusion constants Db,t,p with functions that diverge as g — according to the scaling 
laws 



Db,t (q) = q±' fB,T 



m 



(5.89) 



50 



where we've used the exact d = 2 exponents -2 = | and C = | in the general scahng law 
( ^.65| ). -D||,p, on the other hand, are, like A, constants, since z = 2( (see general equation 
( p.l2| )), which also makes sense since D\\^p are unrenormalized graphically. Hence, the only 
replacement needed to turn the harmonic results into the correct results for the full, non- 
linear theory in ci = 2 is ( [5.89| ). 



D. d > 2 

Now we turn to li > 2. Here the canonical exponents need not hold if 7^ 0. The 
obvious thing to do, therefore, is to determine whether a small A2 is a relevant or irrelevant 
perturbation at the A2 = fixed point. We have attempted to do this to leading (one- loop) 
order in a 4 — e expansion. This involves calculating the perturbative corrections Gi 2 and 
Go-a in the recursion relations ( |5.15|) - ( |5.13| ) to one loop order, and to linear order in A2. 



Once this is done, we can find a fixed point with A2 = 0, and then calculate the linear 
renormalization group eigenvalue of A2 at this fixed point. That is, we'll expand the right 
hand side of the recursion relation ( [5.15| ) for A2 to linear order in A2, obtaining: 



'-^^rA. (5.90) 

for Ai = A*, Xp = A*, (72 = (y\, and A2 small. If 72 is < 0, then the A2 = fixed point is 



(at least locally) stable, and the canonical exponents (|5.80|) - ( ^.82|) will hold for all d. Un- 



fortunately, an (extremely laborious!) calculation (involving 14 different Feynmann graphs) 
shows, after many seemingly miraculous and unexpected cancellations between different 
graphs, that and G0-2 are exactly zero to one loop order. This implies that 

;^ + 1 - ;2 = 0(6^) (5.91) 

and that, to this order at least, Ai^2,p and can take on any value at the fixed point. That 
is, to this order, there appears to be a fixed "line" (actually, a fixed 4 dimensional subspace 
(Ai, A2, Ap, (72)), instead of a single fixed point. This, unfortunately, eliminates all of our 
predictive power for the exponents. For example, keeping D\\ fixed leads to 
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z-2C = -G||(Ai,A2,Ap,a2). (5.92) 

Our earlier arguments show that vanishes if A2 does; however, to this order Ag can be 
anything; hence, so can Gy, and so we get no information about z and C from this relation 
at all. Likewise the recursion relation for A leads to 

z-C-2x+l-c? = -Ga(A2, Ai, Ap, (T2) (5.93) 

with the right hand side again vanishing if A2 = 0, but taking on any value you like if A2 
can be anything, as it can, to this order. 

So what actually happens for 2 < < 4? Unfortunately, from our one loop calculation, 
we cannot say, but can only enumerate the possibilities: 

Possibility I: At higher order, A2 proves to be irrelevant, and flows to zero at the fixed 
point. In this case, the canonical exponents ( [5.80| ) - ( [5.82| ) will hold, for all flocks, in 
all d in the range 2 < d < A. 

Possibility II: ^ = to all orders (i.e., exactly). In this case, there is a fixed line 
(or, more generally, D-dimensional subspace with D > 1) with exponents that vary 
as continuous functions of A2, which can take on any value. Hence, the exponents 
X, and C, will be continuously variable functions of the parameters in the ordered 
phase. This behavior is somewhat reminiscent of that of the d = 2 equilibrium X — Y 
model, although here it is occurring for an entire range of spatial dimension 2 < d < 4, 
and, furthermore, is not associated in any way with the absence of true long ranged 
orientational order, since such order is actually present in our model. 

Possibility III: ^ > at higher order, and the ordered phase is controlled by a new, 
A2 7^ fixed point. In this case, the exponents will again be universal, but presumably, 
different from the canonical ones ( |5.80| ) - ( |5.82| ) . Unfortunately, we have no idea what 
they will be. 
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We should emphasize that all flocks will be described by only one of the three possibilities 
enumerated above (i.e., one can't have different possibilities realized in different flocks). 
Unfortunately, we have no idea which of the above possibilities is realized for 2 < d < 4. 

VI. ANISOTROPIC MODEL 

Not all flocks, of course, are equally likely to move in any direction in the space they 
occupy. Flocks of birds, for instance, although they occupy a. d = 3 dimensional volume (the 
air), are far more likely to move horizontally than vertically. This is presumably because 
gravity breaks the rotational symmetry between the horizontal plane and vertical directions. 

One can imagine a variety of "microscopic" rules, like the Vicsek rule described earlier, 
that would exhibit such anisotropy. For example, one could apply a "Vicsek" rule in three 
dimensions, selecting thereby a vector n. Instead of moving along that vector, however, one 
could instead move along a vector "compressed" along some {z) axis: 

n' = sUzZ + n±_ (6.1) 

with s < 1 and = n — rizZ. This will tend to promote motion in the x — y plane at the 
expense of motion in the 2;— direction. Alternatively, one could project all velocities into the 
x — y plane, apply a Vicsek rule to them (while still sampling neighbors in three dimensions), 
and then add to this xy move a random decorrelated step in the z direction [p!0| . 

For technical reasons that will, we hope, become obvious, we'll focus our attention on 
systems which, whatever their spatial dimension d, have an easy plane of motion; i.e., two 
components of velocity that are intrinsically favored over the other d — 2. We'll also assume 
perfect isotropy within this plane and within the d — 2 dimensional "hard" subspace. The 
case of birds flying horizontally corresponds to d = 3. 

A natural extension of our fully isotropic model (EOM) to this case is 

dtv + Ai {v-V) v + A2 (V ■ v) v + A3V (jvf) = 
-VP(p) + av- /3 \v\'^v- 5avH + DbV (v ■ t/) + D^Vlv 
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+D^Vlv + D2{v-Vyv + f (6.2) 

Mass conservation, of course, still applies: 

dtp + W ■ (pv) = (6.3) 
and the pressure -P(p) will still be given by the same expansion in 6p = p — po 

oo 

P{p) = T.^n{Spr ■ (6.4) 



n=l 



In equation ( p.2| ), vh denotes the d — 2 "hard" components of v, i.e., those orthogonal to 
the d = 2 easy plane. Likewise, Vg and denote the operators Y.'i=i ^ and Y.i=3 
respectively, where i = 1,2 are the "easy" Cartesian directions, and i = 3 ^ d the "hard" 

I — * 1 2 

ones. The term —6a \ vh\ , 6a > suppresses these components relative to those in the easy 
plane. 

Equation (|6.2| ) is not, of course, the most general anisotropic model we could write 



down. For instance, one could have anisotropy in the non-linear terms: e.g, terms like 
(ve ■ Ve could havc different coefficients than (jJh ■ vh- However, because vh winds 
up being "massive," in the sense of decaying to zero too rapidly (ie, non-hydrodynamically) 
at long wavelengths and times to non-linearly affect the hydrodynamic (long wavelength, 
long time) behavior of the flock (in its low temperature phase), any additional terms in 
( |6.2[ ) distinguishing vh and Vg will have no effect on the hydrodynamic behavior in the low 
"temperature" phase. That is, ( |6.2| ) already contains enough anisotropy to generate all 
possible relevant, symmetry allowed terms in the broken symmetry state. Hence, we will 
keep things simple and not generalize ( |6.2| ) further. 

As we did for the isotropic problem, we will now break the symmetry of this model, i.e., 
look for solutions to the form: 

v{f,t) = {v) + 6v{f,t) (6.5) 

Now, however, the direction of the mean velocity (v) (which we'll chose as before, to be 
a static, spatially uniform solution of the noiseless (^f = 0^ version of ( |6.2| ) ) is not arbitrary, 
but must lie in the easy (1, 2) plane. To see this, let's, without loss of generality, write 
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(v) = VoyV + VozZ (6.6) 

with Voy and Voz constants, y in the easy plane and z one of the d — 2 "hard" directions. To 
solve (|6.2| ) with / = 0, these must obey 

aVoy - P {Vly + Vl^) Voy = (6.7) 

and 

{a - 6a) - P « + vl) Voz = (6.8) 

Subtracting Voyy< ( |6.8|) from Wq^x (|6.7|) we obtain 6aVoyVoz = 0, which imphes that either 
Voy or Voz must be zero. It is straightforward to show that the former solution is unstable 
(with two linear eigenvalues a > 0) to small Ve fluctuations, while the latter is stable (with 
d — 2 linear eigenvalues —6a < 0) to vh fluctuations, so the solution with (v) in the easy 
plane is the stable one. Furthermore, fluctuations in the "hard" directions are "massive," 
in the sense of decaying rapidly to zero even at long wavelengths, and so can be neglected 
in the low temperature phase (just like v\\ fluctuations in the isotropic case). Likewise, if we 
take 

{v) = VoV (6.9) 

fluctuations in 6vy = Vy — Vo, will also be massive (with linear eigenvalue— 2q;) . Eliminating 
the massive fields 6vy and vh in favor of the pressure, as we did for 6v\\ in the isotropic case, 
gives 

6Vy = -DpydyP (6.IO) 



vh = -DphVhP , (6.11) 
where we've defined the diffusion constants 

fl„ - ^ (6.12) 
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DpH = ^ (6.13) 

da 

and we've used the relation ( |6.4| ) for the pressure, and dropped all but the leading order 
linear terms in 6p, since higher powers of 6p in Eqns. ( |6.10| ) and ( |6.11| ) prove to be irrelevant. 



Using the solutions (|6.10|) and ( |6.11|) , and taking, for the reasons just discussed. 



V (f) = {vo + 5vy (f, t)) y + v^ (f, t)x + VH (r, t) (6.14) 
we can write a closed system of equations for {f, t) and 6p (r, t): 

dt6p + Vodydp + (pv^) = (Opydl + DpHV]j) 6p (6.15) 



+ + D^dl + DhVI) Vx + U (6.16) 



where we've defined A = Ai + A2, and 7 = Ai^o, and dropped irrelevant terms. 

Proceeding as we did in the isotropic model, we begin by linearizing these equations, 
Fourier transforming them, and determining their mode structure. 

The result of the first two steps is the Fourier transformed equations of motion 

[-i {u - VoQy) + Fp {q)] 6p (g, uj) + iqxPoVx (<?, uj) = (6.17) 

[-i {uj - 'jqy) + T,, (q)] (g, uj) + icTiqxSp {q, uj) = fx (g, uj) (6.18) 
where we've defined 

T,iq} = D,yql + D,Hql (6.19) 

r,{q)=D\^ql + DHq% + Dxql . (6.20) 

Again as in the isotropic model, we first determine the eigenfrequencies uj {q) of these 
equations, finding 
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^± (q) = c± (t)q) q - ie± (q) (6.21) 

where the sound speeds 

c± (^g, 0g) = ^ (7 + ^0) COS % ± C2 (% 0,-) (6.22) 

with 



C2 (% 0g) = y ^ (7 - ^o)^cos2%+ criposin^6'g'cos2 0g , (6.23) 

where 9^ is the polar angle between q and the y-axis, and 0q-is the azimuthal angle, measured 
relative to the x-axis: i. e., the angle between the projection of q orthogonal to y, and the 

X~dbX.\S. 

A polar plot of this sound speed versus for (p^= (i. e., g in the "easy" (i. e., x — y) 
plane) looks exactly like that for the isotropic model (figure 2). Indeed, any slice with fixed 
0q- looks qualitatively like that figure, although, as 0^- — ^ (i. e., as gl, the projection of q 
orthogonal to y, approaches orthogonality to the x-axis), the sound velocity profile becomes 
two circles with their centers on the y axis and both circles passing through the origin. 



The dampings e± (g) in (|6.21|) are (g^), and given by 



e± (g1= ±|^|^(r.(g1 + r,(g-)) 

2C2 {6ci, (l)q) 

T|^^(r.® + ^r,ffl) (6.24) 

2C2 (fc',-, 0g-) V Vo J 

Note that, unlike the isotropic problem in c? > 2, here there are no transverse modes in 
any d: we always have just two longitudinal Goldstone modes associated with 5p and v^- 

We can now again parallel our treatment of the isotropic model and calculate the corre- 
lation functions and propagators. The calculation is so similar that we will not repeat the 
details, but merely quote the results: 

G„„to-,^) = ''"-"°^'>''''^"^ , (6.25) 

Den (g, uo) 



Den (g, u) 



G.,{q,u;)= , (6.26) 
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(6.27) 



A {co-Voqyf + rliq) 



and 



I Den (q, u) 



Cpv [q, a;) = [Sp [q, a;) {-q, -u)) = — J ^ , (6.30) 

\lJen [q, uj)\ 



Cpp{q,u;)= J^ff^ , (6.31) 
\JJen [q, uj)\ 



where we've defined 



Den {q, uj) = {u - c+ (%, (p^) q) (a; - c_ (%, (p^) q) 

+ i [u; (r, (q) + (q)) - qy {voV, (q) + (q))] (6.32) 

which, of course, imphes 

\Den (q, u;)f = (u; - c+ (%, 0^) qf (u - c_ (%, 0^-) g)^ 

+ [a; (Fp (g) + (g)) - (^;,r, (g) + 7!^ (g))]' (6.33) 

These horrific expressions actually look quite simple when plotted as a function of u at fixed 
g; indeed, such a plot of looks precisely like the solid line in figure 6: two asymmetrical 
peaks, centered at a; = c± {9^, (f)^) g, with widths e± (q) oc g^. 

Note that, at this linear order, everything scales as it did in the isotropic problem: peak 
positions oc g, widths oc g^, and heights oc ^. 

Continuing to blindly follow the path we trod for the isotropic problem, we can calculate 
the equal-time — correlation function: 
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Cvv (q) = {vx {q, t) Vx {-q, t)) 

= I '^C^^{q,uj) 
J-oo /vr 

_ A <P{q) 

2 Tl {q) 

where (g) depends only on the direction q of g, and is given by 



(6.34) 



iq) ^ 



C2 (t^g, <Pq) q 



(c+ <pq) q - VoqyY 



c. 



+ K^q, VqJ q - Voqy + (c+ ^ - XlVoqy) ^ 

2 



^ (c_ (^g-, 0,-) g - Voqy) 



c+ 05-) g - Voqy + (c_ (6'g-, (p^ecq) q - y^iVoqy) Y7 



(6.35) 



These fluctuations again diverge hke ^ as |g| — 0, just as in the isotropic problem. 

This completes our abbreviated discussion of the linearized theory of the anisotropic 
model. The most succinct summary of this linearized theory is that everything scales just 
as it did in the isotropic problem. This implies that the non-linearities (i. e., the A and 
(72 terms in the equations of motion (|6.15| )) become relevant in and below the same upper 



critical dimension due = 4 as in the isotropic problem. For c? < 4, therefore, these non- 
linearities will change the long-distance behavior of the anisotropic model. We will now 
treat these non-linearities using renormalization group arguments similar to those we used 
for the isotropic model in d = 2. Now, however, they will work for all d between 2 and 4. 
Notice that all of the non-linearities in ( |6.15| ) are total x-derivatives, just as in the d = 2 



case for the isotropic problem. Now, however, this is true in a// spatial dimensions, not just in 
d = 2. (This, of course, is the reason we chose to consider precisely two "soft" components). 
Thus, we will now be able to derive exact exponents in this model for all spatial dimensions. 
We will not go through the arguments in detail, as they are virtually identical to those in 
the d = 2 case for the isotropic model, but will simply quote the conclusions: 

1. There are no graphical corrections to any of the diffusion constants in ( |6.15| ) except 

2. The stable fixed point that controls the ordered phase must have A* 7^ at least for 
A(0) < Ap(0), which is a finite fraction of all flocks, and 
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3. A and Xp are not graphically renormalized. 



Point one suggests that, in constructing our dynamical renormalization group for ( |6.15D 



we should scale the x-direction differently from both the ^/-direction and the d — 2 hard 
directions. Furthermore, since both the y-direction and the d — 2 hard directions are alike in 
having their associated diffusion constants unrenormalized, we should scale these directions 
the same way. Therefore, in our renormalization group, we will rescale as follows: x bx, 
{y, xh) b'^iy, xh), t — > ¥t. With these rescalings, the recursion relations for Di, i x, p, 
A, and Xp become: 

^ = {z-2C)D, , (t^x), (6.36) 
^ = [z-2x + {l-d)C- 1] A (6.37) 

^ = {x + z-l)Xp (6.38) 

All three relations are exact, since none of these parameters experiences any graphical renor- 
malization. As in the isotropic case, we want all of these parameters to flow to fixed points; 
this leads to three exact scaling relations between the three exponents z, and (: 

z = 2C (6.39) 

z-2x+{l-d)C = l (6.40) 

X = 1 - z, (6.41) 
hose solution is easily found in all d < 4: 
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1-d , , 

Note that these reduce to our isotropic results in (i = 2, as they should, since the two models 
are identical there. They also reduce to the harmonic values z = 2, ( = 1, and x = in 
(i = 4, as they should, since 4 is the upper critical dimension. 
In the physically interesting case of (i = 3, we obtain: 

C = I (6.45) 
z = ^ (6.46) 

X = ~ (6.47) 

As in the isotropic case, we can use scaling arguments here to show that the effect of the 
non-linearities can be fully incorporated by simply replacing everywhere it appears in 
the linearized expressions by the divergent, wavevector dependent scaling form: 



(q) = qlr'f 



Hv.] i M. 
A V A 



(6.48) 



,(f if) 

Doing this leads to all of the scaling laws for this anisotropic problem quoted in the intro- 
duction. 



VII. TESTING THE THEORY IN SIMULATIONS AND EXPERIMENTS 

In this section, we discuss how our theory can be tested in simulations and direct ob- 
servations of real flocks. The "real" flocks may include, e.g., mechanical, self-propelled "go 



carts" packed so densely that they align with their neighbors |]TT|, as well as aggregates of 
genuinely living organisms. 

We begin with a few suggestions about the best boundary conditions and parameter 
values for simulations or experiments, and then describe how the correlation functions and 
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scaling exponents x, z and C, predicted by our theory can be measured. The most useful 
boundary conditions are "torus" conditions; that is, reflecting walls m. d — \ directions, 
and periodic boundary conditions in the remaining direction, call it y (see figure 9). The 
advantage of these conditions is that one knows a priori that, if the flock does spontaneously 
order, its mean velocity will necessarily be in the periodic (?/) direction. 

It might be objected that imposing such anisotropic boundary conditions breaks the 
rotation invariance our model requires, but this is not, in fact, the case. A "bird" deep 
inside the box moves with no special direction picked out a priori; it can only flnd out about 
the breaking of rotation invariance on the boundary if the bulk of the flock spontaneously 
develops long range order. This is precisely analogous to the way one speaks of a ferromagnet 
as spontaneously breaking a continuous symmetry even if it orders in the presence of ordered 
boundary conditions. 

So, by imposing these boundary conditions, we know the direction of the flock motion 
(the y direction in the simulation), and, therefore, have oriented the simulation axes with 
the axes used in our theoretical discussion; i.e., our || axis equals the simulation's periodic 
direction. 

Alternative boundary conditions add the additional complication of having to flrst de- 
termine the direction of mean flock motion before calculating correlation functions. This 
complication is even worse for a flnite flock (as any simulation must treat), since the mean 
direction of motion will wander, executing essentially a random walk that will explore the 
full circle in a time of order Tji^ck = ^tt^^. Our results, which assume a constant direction 
of flock motion, will only apply for time scales t « Tpiock- Even drifts of the mean flock 
direction through angles << 27r can cause problems, however, since most of the interesting 
scaling behavior is concentrated in a narrow window of angles gy ^ q]_ ^ q±', i.e., near the 
direction of mean flock motion. So this drift greatly complicates the experimental analysis, 
and is best avoided by using the toroidal boundary conditions just described. 

Of course, it is considerably harder to produce these boundary conditions in a real 
experiment. Ants walking around a cylinder may come close, although gravity will always 
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break rotation invariance on a real cylinder. Perhaps the experiment could be done on the 
space shuttle, or with a rapidly spinning cylinder producing artificial gravity that swamps 
real gravity, or by using neutrally buoyant organisms in a fluid. Alternatively, one could use 
a "track" such as that shown in figure 10, and take data only from the cross-hatched region, 
chosen to be in the middle of the straight section of the track, far from the curves. 

Other, more ingenious ways to pre-pick the direction of mean motion through boundary 
conditions may also be dreamed up by experimentalists more clever than we are. 

We strongly caution anyone attempting to test our results, however, that it is only 
through boundary conditions that one may prepick the direction of mean motion. Any 
approach that prepicks this direction in the bulk of the flock, such as giving each bird a 
compass, letting them be blown by a wind, or run downhill, or follow a chemical scent, 
etc., will lead to a model outside the universality class of our isotropic model, since the 
starting model does not have any rotation invariance to be spontaneously broken (unless 
the anisotropy leaves an "easy plane" in which all directions are equivalent, in which case our 
anisotropic model of section VI applies). Indeed, such flocks of "birds with compasses" will 
be less interesting than the models we've studied here, since the "compass" will introduce a 
"mass" that makes any fluctuation away from the pre-picked direction of flock motion decay 
rapidly (i.e., non-hydrodynamically) with time. In such a model, it's easy to show that the 
non-linearities are irrelevant, and there are no interesting fluctuations left at long distances 
and times. 

And now a few words about parameter choices. For definiteness, we will discuss in what 
follows the Vicsek model, whose parameters are Vq — where S is the distance the birds 
travel on each time step and Rq is the radius of the circle of neighbors, the mean number 
density po in units of where d is the dimension of the system, and the noise strength 
A, which is the mean squared angular error. Since the interesting non-linear effects in our 
model come from terms proportional to Vq, those effects will become important at shorter 
length scales in a faster moving flock. That is, in, e.g., the Vicsek model, should we chose the 
dimensionless velocity as large as possible, consistent with the flock ordering. However, if we 
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take vq too big, i.e., ^ 1, then, on each time step, each bird is hkely to have a completely 
different set of neighbors. It is difficult to see how order can develop in such a model. So, to 
take vq as big as possible without violating vq ^ 1, we should chose fo ~ 1. The simulations 
of Vicsek et al. |I| took fo ^ 1, and, hence, probably never explored (in their finite flocks) 
the long length scale regime in which our non-linear effects become important. 

Now to the mean density po, which is, of course, just determined by the total number 
of birds and the volume V of the box via po = y. We clearly want this to be large 
enough that each bird usually finds some neighbors in its neighbor sphere: This means we 
want PqRq > 0(1). However, if we make po too large, each bird has so many neighbors 
that a simulation is considerable slowed down, since the "direction picking" step of the 
Vicsek algorithm takes a time proportional to the number of neighbors (because we've got 
to average their directions). Thus, for simulations, one wishes to choose po as small as 
possible, consistent, again, with getting good order. 

Finally, we consider the noise A. Here again, to see our fluctuation effects, we want A 
as big as possible. However, if A is too big, the flock won't order. Furthermore, even if A 
is small enough that the flock it does order, we want also to be sure that we are well below 
the critical value Ac of A at which the flock disorders. Otherwise, for distances smaller than 
the correlation length ^ associated with the order-disorder transition, the scaling properties 
of the flock will be controlled by the fixed point that controls the order- disorder transition, 
not the low temperature fixed point we have studied here. 

If this transition is continuous, as it appears to be in Vicsek's simulations this cor- 
relation length diverges as A ^ A^. Thus, to observe scaling behavior we predict over as 
many decades of length scale as possible, we want to choose A substantially less than Ac, 
but as big as possible consistent with this (to maximize fluctuation effects). Choosing A 
to be a little below the point at which the mean velocity (v) starts to "saturate" seems 
like a fairly good compromise between these two competing effects. Similar considerations 
apply for choosing the optimal po, and vq, which we want to be as small or big, respectively, 
as they can be without substantially suppressing long ranged order. The best choices will 
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probably lead to all three parameters po, "^o and A begin, in suitably dimensionless units, 
0(1). 

Having chosen the appropriate parameter values and boundary conditions, what should 
an experimentalist or simulator measure to test our theory? 

We have already discussed a number of such measurements in the introduction; namely, 
the spatially Fourier transformed equal-time and spatio-temporally Fourier transformed un- 
equal time density-density correlation functions Cp{q) and Cp{q,uj), respectively. Our 
predictions for these are given in Eqns. ( |5.38|) and (|5.4CI| ). 



One additional correlation function that can be measured quite easily is the mean squared 
lateral displacement of a bird 

w\t)^(j4{t)-xm\"'j (7.1) 

perpendicular to the mean direction of motion of the flock. This can easily be measured as a 
function of time in a simulation or experiment simply by labeling a set of n birds in a "strip" 
near the center of the channel with its long axis running parallel to the mean direction of 
bird motion (see fig. 11) and then following their subsequent motion. It is best to center 
the strip in the channel so as to postpone the birds reaching the reflecting walls as long as 
possible. Once they do reach the walls, of course w'^{t oo) saturates at ~ L\, L± being 
the width of the channel. We will deal in the following discussion with times much smaller 
than that required for a bird at the center of the channel to wander out to its edge. Since 
the mean x±— position xf- of each bird obeys 

4{t)=xi{0)+ fvt{t)dt (7.2) 

where v^{t) is the _L velocity of the i'th bird at time t, the mean width is given by 

w\t) = f dt' f dt" {vtit') . vtit")) . (7.3) 



Now we need to relate the velocity of the i'th bird to the position and time dependent 
velocity field v_[_ {r,t). This is easily done: 
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v^{t) = v^m),t) (7.4) 
where fi{t) is the position of the i'th bird at time t. This is given by 

= fi{0) + vtx\\ + 5xi{t)x\\ + 5x^{t) (7.5) 

where 



(7.6) 



is the velocity averaged over all birds, which, as discussed earher, is not to be confused with 
the space averaged f o = j v (r, t) d'^r that appears in the expression for the sound speeds 
c± (%)• This distinction proves to be crucial here, as we shall see in a moment. In ( [7.5[ ), 
5x\{t) and 5x^()f:) reflect the motion of the i'th individual bird relative to the mean motion 
of the flock (at speed v). 

Using (|7.5|) , we see that the desired single bird autocorrelation function in (|773|) is 

{vHt') ■ {t")) = 

{v± (fi(0) + {vt' + 6x11 it')) h + (f) , tf) ■ (fi(0) + vt" + 6x\i (t") X|| + 6x± {t") , t")) 

= Cc {x± {t') - x± {t") , v\t' - t"\ + 6x11 it') - ^^W it") . t' - t") (7.7) 

where Cc (^R, tj is the real space velocity field auto-correlation function defined in the intro- 
duction. 

We assume (and will verify a posteriori) that both 6x\\ and 6xj_ are small enough com- 
pared to the average motion vtx\i that their effect on the velocity- velocity autocorrelation 
in (|7.7| ) is negligible. For now neglecting them, we see that we are left with the task of 
evaluating (^R_i_ = 0, R\\ = vt, t^ . 

Expressing Cc in terms of its Fourier transform then gives 

Cc {R± = 0, R\\ = vt, t) = J d'^-\^dq\\du e<^-'^^ii*)ai (g, u) (7.8) 

Using the fact that Cu (g, u) is peaked a.t u = c± {6^) q with widths that scale like 
((^^''iV)' dominating peak is at u; = uj- for fo(0) > 7(0) or at cij = u;+ for 
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vo{0) < 7(0), with heights that scale hke q^^g (j^j^^ with S = 2x + z + C + d- l 
(see Eqn. ( [5.56| )). Assuming that fo(0) > 7(0), it is straightforward to show that, upon 
integrating ( |778| ) over u, we obtain 

a = / d'^-ig^cig||e<^-K>-^^ii)V- (j^) qt' (7.9) 

This integral is dominated, as t — >• 00, by q\\ ~ (gj.^o)'' /^o ^ (lu hence, 0, and we 

get 

= / rf'^-igx^^gye^^^^-'^'^iiV- (t^^I • (7.10) 
We can scale the time dependence out of this integral with the change of variables 

gi, = ^ and gl = ^ (7.11) 

which give 

C, oc t^^^/^ . (7.12) 

Using this in ( |7.7D for the single bird velocity autocorrelation function, and then using that 
autocorrelation function in the expression ( [7.3|) for the mean squared random walk distance 
gives 

w\t) oc r dt' t dt" \t' - fl^""'^ (7.13) 

JO JO 

Now, we need to distinguish 2 cases: 

Case (1): ^ > — 1 . In this case, which holds in (i = 2, where x = ~| and C = |, the 
double integral over t' and t" is dominated, for t ^ to^ the microscopic time scale, 
by t', t", and \t' — 1"\ of order t ^ tg- Hence, our calculation of Cc which used the 
hydrodynamic (i.e., long time) limiting forms of the correlation functions is correct, 
and ( [7.13| ) holds. Changing variables to T' = j and T" = y, we see that 

w\t)oct'M) =t'/^ ^'^■^^^ 
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the last equality holding in d = 2. Note that this behavior is "hyperdiffusive" : the 
mean squared displacement w'^{t) grows faster than it would in a simple random walk; 
i. e., /aster than linearly with time t. 

Case (2): ^ < — 1 . In this case, which certainly holds for c/ > 4 (where x = 1 ~ f < ~1 
and ( = 1), the integral over t" converges as \t' — 1"\ oo. Hence, that integral is, in 
fact, dominated by \t' — t"\ = 0{to), the microscopic time, where our hydrodynamic 
result Eqn. ( [7. 121 ) is not valid. Presumably, the correct \t' — t"\ — > limit of the single 
bird velocity autocorrelation ( |7.3D is finite; and, hence, so the integral over t" in ( [7. 13[ ) 
approaches a finite limit as t — oo. 

Hence, we get 

w (t) (X / dt' X finite constant oc t , — < — 1 . (7-15) 
Jo C 

We now need only verify our a posteriori assumptions that 5x\\ and 5x± were negligible 
in the velocity- velocity autocorrelation Eqn. ( [7.7| ). 

First consider x±; we have just shown that the root-mean-squared \x± {t') — xj_ {t")\ oc 
\t' — t"\^^~. iFiom our scaling expression (6.43), we see that Cc (^R±, R\\,tj ^ 
= 0,i?||,t) ifi?i<i?||. In (|73), we are interested in _Rj_ oc \t' — anc 



R\\ oc \t' — hence, the condition <^ R\\ will be satisfied as \t' — 1"\ oo pro- 
vided C + X < 1- Since C ^ 1 ^ind x < for all d > 2, this condition is satisfied for all 
d > 2. For 6x\\ we need only show that 6x\\ {t') — fey {t") ^ \t' — t"\ as \t' — t"\ oo. 
This is easily shown by using the fact, alluded to earlier, that 6v\\ (r, t), the fluctuation 
of the velocity along the mean direction of motion, has only short ranged temporal 
correlations. Using this fact, it is straightforward to show that Sx\\{t) just executes a 
simple random walk; that is 



(f) - 6x\\ (t")f oc y^(t' - t") < \t' - t"\ (7.16) 



and hence these fluctuations are negligible as well. 
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Unfortunately, the analogous calculation for the anisotropic model shows that this 
random "transverse walk" is much less interesting: the mean squared transverse dis- 
placement in the x-direction (the direction in the "easy place" of the anisotropic model 
orthogonal to the mean direction of motion, y) is given by an expression very similar 
to ( TO 



(|x,(t) - x,(0)|') = w\t) = [dt'jyt" {vr. it')v,, it")) (7.17) 



and a calculation so closely analogous to that just given for the isotropic model that 
we shan't bother to repeat it for this case shows that 

{v,,{t')v,,{t"))oc\t'-t'f-''--^ (7.18) 

As in our analysis of the isotropic case, here, too, the question of whether "simple 
random walk" behavior (u'^(t) oc t) or "hyper diffusive" behavior (u'^(t) (xP,'y > 1) 
occurs hinges entirely on whether the exponent in ( |7.18| ) is greater or less than —1, 
with hyperdiffusive behavior occurring in the former case (exponent > — 1) and simple 
random walk behavior in the latter (exponent < —1). Using our exact result ( |6.43| ) 
for z in the anisotropic model for 2 < d < 4, we see that hyperdiffusive behavior will 
occur if 

2 2 — 2d 

3-d-- = ^^>-l (7.19) 
z 3 

which is satisfied only for d < 5/2. Unfortunately, this condition is not satisfied 
for either d = 3 or d = 4. In d = 2, the anisotropic model is the same as the 
isotropic model, while for d > A, z = 2 and 3 — — | < — 1. So in no case in which 
the anisotropic model is different from the isotropic one is hyperdiffusive behavior 
observable; rather, we expect w'^{t) oc t for all those cases. This negative prediction 
could be checked experimentally, although its confirmation, while a non-trivial check 
of our theory, would clearly be less exciting than verification of our hyperdiffusive 
prediction w'^{t) oc t^/^ for the isotropic d = 2 model. 
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Some of the numerical tests discussed in this section have been carried out recently, and 
good agreement with our prediction has been reached. |12| 



VIII. FUTURE DIRECTIONS 

In this paper, we have only scratched the surface of a very deep and rich new subject. We 
have deliberately focussed on the most limited possible question: what are the properties of 
a flock far from its boundaries, and deep within its ordered state? Every move away from 
these restricting simplifications opens up new questions. To name a few that we hope to 
address in the coming millennium: 

1. The transition from the ordered (moving) to disordered (stationary, on average) phase 
of the fiock. This can be studied by analyzing the (unstable) fixed point at which 
the renormalized a of our original model ( p. 61) is zero. The dynamical RG analysis 
of this point would be technically similar to the one we've presented here for the low 
temperature phase, with a few crucial differences: 

(a) All components of w, not just the ± components, become massless at the transi- 
tion. 

(b) The fixed point will be isotropic, since no special directions are picked out by (-u), 
since ({7) still = at the transition. 

(c) The (5\v^ V term becomes another relevant vertex. We know, by power counting, 
that at the transition, this vertex becomes relevant in c? = 4. Indeed, if we 
ignore the A vertices, our model simply reduces to a purely relaxational time- 
dependent Ginsburg Landau (TDGL) model for a spin system with the number 
of components n of the spin equal to the dimension d of the space those spins live 
in. 

We have convinced ourselves by power counting that at the transition, for (i < 4, the 
A vertices are a relevant perturbation to the Gaussian critical point. Whether they 
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constitute a relevant perturbation to the 4 — e TDGL fixed point, thereby changing its 
critical properties, can only be answered by a full-blown dynamical renormalization 
group analysis. 

Obviously, a similar analysis could also be done for the anisotropic model. 

2. The shape and cohesion of an open fiock, and its fiuctuations. We have thus far 
focussed on fiocks in closed or periodic boundary conditions. Real fiocks are usually 
surrounded by open space. How do they stay together under these circumstances? 
What shape does the fiock take? How does this shape fiuctuate, and is it stable? 

This issue is somewhat similar to the problems of the shapes of equilibrium and growing 
crystals (e.g., facetting, dendritic growth). In those problems, it was important to first 
understand ftn/A; processes (e.g., thermal diffusion in the case of dendritic growth) before 
one could address surface questions (e.g., dendritic growth). The non-trivial aspects of 
the few/A; processes in flocks (e.g., anomalous diffusion) will presumably radically alter 
the shapes and their fluctuations. 

3. A somewhat related question is: what happens if birds move at different speeds? By 
"move at different speeds", we don't mean simply that at any instant, different birds 
will be moving at different speeds (a possibility already included in our "soft spin" 
dynamical model equation ( p.6|) ). Rather, we mean a model in which some birds 
have a different probability distribution of speeds than others. (In our model, this 
distribution of the speed of any given bird is the same over a sufficiently long time, 
and controlled by the values of the parameters a and /3, with large a and (5 leading to 
a distribution sharply peaked around a mean speed f o = y^, while small a and /? lead 
to a broader distribution). More generally, one could imagine 2 (or many) different 
species of birds, (labeled by k) all flying together, each with different mean speeds v^. 
What would the hulk dynamics of such a flock be? Would there be large scale spatial 
segregation, with fast birds moving to the front of the flock, and slow birds moving to 
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the back? If so, how would such segregation affect the shape of the flock? Would it 
elongate along the mean direction of motion? Would this elongation eventually split 
the flock into fast and slow moving flocks? 

4. At the other extreme, one could consider flocks in confined geometries; e.g., inside a 
circular reflecting wall in d = 2. In such a case, the time averaged velocity of the flock 
{v{f,t))^ could not be spatially homogeneous but would have to circulate around the 
center of the circle; i. e., {v{r,t))^ = f{r)6. The spatially inhomogeneous pattern of 
velocity and density that resulted could be predicted by our continuum equations. This 
problem is potentially related to the previous one, since one way a flock containing, 
say, some very fast birds and other very slow birds, could stay together would be for 
the fast birds to fly in circles inside the essentially stationary volume of space filled 
by the slow birds. It would be very interesting to make the connection between our 
continuum theory and the recently observed circular motion of Dictyostelium cells in 
a confined geometry. 



5. One could relax the constraint on conservation of bird number, by allowing birds 
to be born, and die, "on the wing". Numerical studies of such models, which may 
be appropriate to bacteria colonies, where reproduction and death are rapid, as well 
as the migration of, e. g., huge herds of caribou over thousands of miles and many 



months, have already been undertaken |]14[; it should be straightforward to modify 



our equations by adding a source term to the bird number conservation equation. 

6. It is possible that phase transitions other than that from the moving to the non-moving 
state occur in flocks. For example, in some preliminary simulations of microscopic mod- 
els in which birds try to avoid getting too close to their neighbors, rather than merely 
following them, we have observed (literally by eye) what appears to be a "flying crys- 
tal" phase of flocks: the birds appear to lock themselves onto the sites of a crystalline 
lattice, which then appears to move coherently. It would be very interesting to test nu- 
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merically whether this optical appearance reflects true long-range translational order, 
by looking for a non-zero expectation value of the translational order parameters. 

p^{t) = (s,e^^->W/iV) (8.1) 

which will become non-zero in the thermodynamic (A^ — > oo) limit at a set of re- 
ciprocal lattice vectors G if such long ranged order actually develops. It will also be 
extremely interesting to include the possibility of such long ranged order in our analytic 
model, and study the interplay between this translational order and the anomalous 
hydrodynamics that we've found here for "fluid" flocks. Will anomalous hydrodynam- 
ics suppress the "Mermin- Wagner" fluctuations of translational order, just as it does 
those of orientational order, and lead to true long ranged translational order, even in 
d = 27 Will the crystallization suppress orientational fluctuations, and thereby slow 
down the anomalous diffusion that we found in the fluid case? And in any case, what 
are the temporal fluctuations of PQ{t)1 

It should be noted that this problem potentially has all the richness of liquid crystal 
physics: in addition to "crystalline" phases, in which the set {G} of reciprocal lattice 
vectors in ( |8.1| ) spans all d dimensions of space, one could imagine "smectic" phases in 
which all the G"s lay in the same direction; and "discotic" phases in = 3, in which 
the G's only spanned a two-dimensional subspace of this three-dimensional space. 
The melting transitions between these phases and the "fluid" , moving flock, as well as 
possible direct transitions between them and the stationary flock phase, and between 
each other, would also be of great interest. 

We should point out here that these models differ considerably from recently considered 



models of moving flux lattices |]T5| and transversely driven charge density waves p5|JT6 
in that here, the direction of motion of the lattice not picked out by an external driving 
force, but, rather, represents a spontaneously broken continuous symmetry. 

7. Finally, we'd like to study the problem of the growth of order in flocks. This is a 
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phenomenon we've all seen every time we walk onto a field full of geese: eventually, 
our approach startles the geese, and they take off en masse. Initially, they fly in random 
directions, but quickly the flock orders, and flies away coherently. The dynamics of 
this process is clearly in many ways similar to, e. g., the growth of ferromagnetic order 
after a rapid quench from an initial high temperature Tj > T^, the Curie temperature, 
to a final temperature Tf < Tc, a problem that has long been studied and proven 



to be very rich and intriguing. In flocks, where, as we've seen, even the dynamics of 
the completely ordered state is very non-trivial, the growth of order seems likely to be 
even richer. 

Even this list of potential future problems, representing, as it does, probably another ten 
years of research for several groups, clearly represents only a narrow selection of the possible 
directions in which this embryonic field can go. We haven't even mentioned, for example, 
the intriguing problem of one-dimensional flocking, with its applications to traffic flow (and 
traffic jams), a topic clearly of interest. This problem has recently been studied [|1^] and 
found to also show a non-trivial phase transition between moving and non-moving states. 

We expect flocking to be a fascinating and fruitful topic of research for biologists, com- 
puter scientists, and both experimental and theoretical physicists (at least these two!) for 
many years to come. 
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FIGURES 

FIG. 1. A snapshot of a simulated flock that has reached a statistically steady state. Note the 

enormous fluctuations in the density. Quantitatively, the statistics of the spatial Fourier transform. 
Cp (q) obtained from this picture agree with our quantitative prediction equation (1.10). 

FIG. 2. Polar plot of the direction-dependent sound speeds c± (0^), with the horizontal axis 
along the direction of mean flock motion. 

FIG. 3. Plot of the damping Imu versus qy = \q±\ where q± is the projection of wavevector 

q perpendicular to the direction of the mean flock velocity {v} for fixed projection qx of q parallel 

z 

to {v). Note that, for small the crossover between /mo; oc qy and Imuj oc occurs only for 
directions of propagation q very nearly parallel to the mean flock velocity (i/), since C < 1- 

FIG. 4. Plot of the spatio-temporally Fourier-transformed density correlation function Cp (g, uS) 
versus a; for fixed q. It shows two sharp asymmetrical peaks at w = c± (6*^) q associated with the 
sound modes of the flock, where c± (9^) are the sound mode speeds. The widths of those peaks are 
the second mode dampings Imuj± {9^) oc qj_f± 

FIG. 5. Geometry of the anisotropic model. Birds prefer to fly in the "easy" x — y plane. 
We take their (spontaneously chosen) direction of motion within that plane to be y. The in-plane 
direction perpendicular to that is x. In general d, there are d—2 "hard" directions fn perpendicular 
to this easy plane. The anisotropy of scaling is between x and the other d — 1 directions y, fn- 

FIG. 6. Plot of Cll (q,^) and Ctt (q,^) versus u for identical fixed q. Note the smallness of 
the overlap between the transverse and longitudinal peaks. 

FIG. 7. Feynmann graph renormalizing the noise correlations when A2 = 0. There is a factor 
of the external momentum \q±\ = q^ associated with each vertex; hence this graph does not 
renormalize A, but, rather, only changes 0{q^) pieces of the / — / correlation function. 



77 



FIG. 8. Feynmann graph for diffusion constants. When A2 = 0, this graph is proportional to 
at least one power of \q±\ = Qx, and, so, cannot renormalize Dy or Dp. 

FIG. 9. Illustration of the optimal boundary conditions for simulations and experiments to test 
our predictions. The top and bottom walls are reflecting, while periodic boundary conditions apply 
at the left and right walls (i.e., a bird that flies out to the right instantly reappears at the same 
height on the left). The mean direction of spontaneous flock motion, if any occurs, is clearly forced 
to be horizontal by these boundary conditions. In spatial dimensions d > 2, one should choose 
reflecting boundary conditions in cZ — 1 directions, and periodic in the remaining direction, thereby 
forcing (v) to point along that periodic direction. 

FIG. 10. More practical "track" geometry for experiments on real flocks. Data should only be 

taken from the cross-hatched region centered on the middle of the "straightaway". 

FIG. 11. Illustration of the experiment to measure the mean squared lateral wandering w'^{t). 
One labels all of the birds some central stripe (of width <C L, the channel width), and then 
measures the evolution of their mean displacements x±^{t) perpendicular to the mean direction of 
motion (which mean direction is horizontal in this figure). 
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